What is a good distance measure for matching SIFT descriptors depending on the distribution of their noise? I have read some papers about distance measures like Euclidean, Manhattan or Chi-Square for matching gradient based image descriptors like those computed from the SIFT Algorithm (128-D vectors). Most of them state, that one or the other measure is more suitable for that task and that this decision depends on the noise distribution assumed (Like L2 is suitable if Gaussian noise is assumed or L1 for Laplacian noise).

So my question is: how do I determine the distribution of my descriptors and how does that affect the choice for a distance measure?

In one paper the authors stated that the distribution would be the distances between SIFT descriptors of correctly matching keypoints. But doesn't that mean this distribution is biased by the choice of the distance metric itself?!

For obtaining distributions or modes in your, a basic approach is to cluster the data. That will of course bias the distribution towards the metric, or in other words the distribution will be tied to the metric. For example, it is possible to quantize the SIFT descriptors by using vocabulary trees to visual words and then plot the 1D histogram of the individual words. That's only one possible way to do this from a large collection of SIFT descriptors. Another approach would use nonparametric methods (Parzen Window, Meanshift...). In either case the distribution will be data-dependent and quite complex.

I think this question also deserves to mention the Root-Sift paper, where authors propose to use the Hellinger distance instead of Euclidean. The implementation is very easy (just take a sqrt) and the results are consistently improved. K-means can still be used. Apparently, Hellinger is a better metric than Euclidean for SIFT.


Matlab frequency response I'm trying to make a frequency response of that function:

in matlab a made this:

d=[1 -0.494 0.6404];
p=[1 1.58567330897 0.96];
[h,w] = freqz(p,d,40000);
z = abs(h);
plot(w/pi,z);

d and p i calculated.
The plot gives me this:

but in book has that answer:

as you can see the peak of magnitude is 2 and not about 6.5. '-'
What am I doing wrong?
H2 has an abs 1 because it is a all pass filter.

Your answer is correct. Please note that your graph is in the range $0$ to $1$ while the graph in the text is in the output range of $-\pi$ to $+\pi$

There is some output scaling issue. But since you are getting the same curve, you are just messing some constants is all.

please try $\text{plot}(w,z)$


Properties of Up- and Downsampling Given a gaussian white noise with mean 0 and $\sigma = 1$.
I'm interested in the PSD (power spectral density), when the signal is plain up- or downsampled. No interpolation, no filtering nothing.

There are some related topics in the web, but they are always ending with someone considering low- or passband filters. They are also often speaking about the aliases which get e.g. back in the nyquist band when downsampling, but they never properly describe it with math:

So, what happens with the given PSD, (which should have an amplitude of 1/fs between -fs/2 and fs/2) when it's led trough a down and upsampler?

I assume first, up and down or the opposite should result in the same PSD. Otherwise one could just up and downsample to benefit from some noise reduction.
I could also imagine that up and down sampling make no changes at all to white noise. On the other hand, there might be an effect when downsampling, due to aliases.

Let me try to explain what happens to the PSD (power spectral density) of a discrete-time Gaussian white noise $x[n]$ of power $\sigma_x^2$ after it's either expanded by $L$ (insertion of L-1 zeros between consequitive samples of $x[n]$), or compressed by $M$ (aka decimation, selecting every M-th sample of $x[n]$).

First, let's put a convenient framework, by defining required quantities and their assumptions.

A discrete-time random process (RP) is constructed by a set of random variables (RV) and represented by $\{X[n,s)\}$ where integer $n$ denotes time-indexing for each RV and $s$ denotes outcome of a random experiment associated with those RVs or the RP. Customarily, an outcome $s$ selects a complete discrete-time waveform $x[n]$, which is an instance of the RP, from the ensamble set of it. The ensamble set includes all possible discrete-time waveforms. For notational simplicity, the sequence $x[n]$ is referred to as the random process, unless stated otherwise.

Probabilistic analysis of the RP $x[n]$ proceeds by defining various joint PDFs (probability density/distrbution functions) between random variables of the RP. Among those joint PDFs, (and moments associated) two of them become most useful and important for practical engineering applications. They're the mean and the Auto-Correlation Sequence (ACS) of the process, as defined below for a real $x[n]$ :

1- Mean of a RP is: $\mu_x[n] = \mathcal{E}\{x[n]\} $

2- ACS of a RP is: $\phi_{xx}[n,m] = \mathcal{E}\{x[n]x[n+m]\}$

As can be seen, in general the mean and ACS of a RP depend on time index $n$. The ACS depending also on the lag m. A very important simplification for the analysis of RP is made by defining what's known as a wide sense stationary (WSS) RP, whose mean and ACS are defined to be independent of time; i.e.,

1- Mean $\mu_x[n]  = \mu_x = \mathcal{E}\{x[n]\} $

2- ACS $\phi_{xx}[n,m] = \phi_{xx}[m] = \mathcal{E}\{x[n]x[n+m]\}$

Now, for given a WSS random process $x[n]$, its ACS $\phi_{xx}[m]$ holds a very important position. The power spectral density (PSD), $S_{xx}(e^{\omega})$, of a WSS RP is defined to be the discrete-time Fourier transform (DTFT) of its ACS:

$$ S_{xx}(e^{\omega}) = DTFT \{ \phi_{xx}[m] \} = \sum_{m=-\infty}^{\infty} \phi_{xx}[m] e^{-j\omega m} $$

At this point, it's clear that the PSD of a RP exists iff the RP is WSS and further that its ACS is stable so that its DTFT converges. I assume  @DilipSarwate may comment if a (major) ill-statement was presented up to here.

we are now ready to tackle your problem !

The compression by integer M (which you refer to as downsampling without filtering) operation is:

$$ x[n] \longrightarrow \boxed{ M \downarrow} \longrightarrow y[n] ; ~~~y=x[Mn] $$

Given an i.i.d (independent, identically distributed) zero mean,  WSS, (Gaussian) white noise RP $x[n]$ with power $\sigma_x^2$ it can be shown that its ACS and PSD are:

$$\phi_{xx}[m] = \sigma_{x}^2 \delta[m]$$
$$S_{xx}(e^{j\omega}) = \sigma_x^2 ~~~, ~~~\text{ for all } \omega $$

Then what is the PSD associated with the compressed output $y[n]$ ? To show this, first we must find its ACS and see that the RP $y[n]$ is indeed WSS.

$$\phi_{yy}[n,m] = \mathcal{E}\{y[n]y[n+m]\} = \mathcal{E}\{x[Mn]x[M(n+m)]\} = \mathcal{E}\{x[Mn]x[Mn+ Mm)]\} = \phi_{xx}[Mm] = \sigma_x^2 \delta[Mm] = \sigma_x^2\delta[m]$$

Hence we see (with a slightly bold step of treating $Mn$ the same as $n$ alone, which's justified by the WSS and whiteness of the input $x[n]$) that the ACS associated with the compressed signal $y[n]$ is identical to the ACS of the input $x[n]$. Hence their PSD's are also the same:

$$S_{yy}(e^{j\omega}) = \sigma_x^2 ~~~,~~~\text{ for all } \omega  $$

Next we apply the same for the expansion (you refer upsampling) process which can be shown to be as:
$$ x[n] \longrightarrow \boxed{ L \uparrow} \longrightarrow y[n] = \begin{cases} { x[n/L] ,~~~n=rL \\ 0 ,\text{ otherwise} }\end{cases} $$

Let's try to compute the ACS of the expanded output $y[n]$.
$$\phi_{yy}[n,m] = \mathcal{E}\{y[n]y[n+m]\} = \begin{cases} { \sigma_x^2 ,~~~ m=0 ~~~\text{ and } n=rL \\ 0 ,~~  m \neq 0~~~\text{ or }  n \neq rL }\end{cases} $$

When computing $\phi_{yy}[n,m]$, I assumed that the stuffed zeros were deterministic random variables with a pdf of $f_X(x)=\delta(x)$, which can only take a value of $0$ with a probability of one. They are not independent in themselves but independent from nonzero samples. Furthermore, since their value is always zero, when multiplied inside an expectation operator they always yield a zero result.

Now there is a problem! Eventhough the input $x[n]$ is WSS, the output is not. As its ACS $\phi_{yy}[n,m]$ depends not only on the lag $m$, but also on the exact time instance $n$. Hence its ACS is not of the form that's suitable for the DTFT expression. So at this point we can say that the PSD associated  with $y[n]$ does not exist in the form of $$S_{yy}(e^{j \omega}) = DTFT\{ \phi_{yy}[m] \} $$

May be a parametric Fourier transform or a Short-Time Fourier transform (STFT) or a 2D Fourier transform can be utilized but I'm not very sure about it.

Finally, we can state the following observation that since the expansion and compression operators are not LTI systems (they are time varying), then we are not surprised to see that their outputs are not WSS even if their inputs are so, which would not be possible with LTI systems.


Difference equation and Impulse response of an IIR filter - deconvolution Given the input to the IIR filter

$$
x[n]=\delta [n]+0.25\delta[n-2205], \quad \text{where}\quad \delta[n]= \begin{cases}1, &n=0\\0, & n\neq 0\end{cases}
$$

What I want to do is to design an IIR-filter that cancels the "echo", such that I'm only left with the $\delta[n]$. To do so I of course make an IIR-filter, where you take the $\mathcal Z$-transform of the impulse response, and cascade it with the inverted

$\mathcal Z$-transform of the input:

$$
X(z)=1+0.25z^{-2205}
$$

Transfer function of the IIR filter:

$$
H(z)=\frac{1}{1+0.25z^{-2205}}
$$

Then our output naturally is:

\begin{align}
Y(z)&=\left(1+0.25z^{-2205}\right)\cdot\frac{1}{1+0.25z^{-2205}}\\
& = 1\\
\Rightarrow y[n]&=\delta [n]
\end{align}

So the transfer function of my filter is $H(z)$. What I want to do is to find the difference equation for my filter, and possibly find the mathematical form of the impulse response. Since the order is so high I'm having difficulties.

I'm pretty sure the difference equation will not be infinite, while the impulse response will... why??

So your transfer function is

$$H(z)=\frac{1}{1+az^{-L}}\tag{1}$$

with $a=0.25$ and $L=2205$. If $Y(z)$ and $X(z)$ are the $\mathcal{Z}$-transforms of the output and input signal of the filter, respectively, then you have

$$Y(z)=H(z)X(z)$$

which, using (1), is equivalent to

$$Y(z)(1+az^{-L})=X(z)\tag{2}$$

Transforming (2) back to the time domain results in the difference equation of the recursive filter:

$$y[n]+ay[n-L]=x[n]\tag{3}$$

You can get the corresponding impulse response either from the inverse $\mathcal{Z}$-transform of (1), or from setting $x[n]=\delta[n]$ in (3) and observing what happens to the output. In either case, the result should be

$$h[n]=\begin{cases}(-a)^{n/L},&& n=0,L,2L,3L,\ldots\\
0,&& \text{otherwise}\end{cases}$$


Roberts Edge Detector how to use? I am trying to use Roberts edge detection to process an image. Do I just apply both of the masks to the image and perform convolution as normal? Could someone give me the breakdown of how to use this edge detection method, as I am trying to program it to process a greyscale image. I convoluted the image using both kernels separately  but the image dent' look right.

Thanks.

Robert's Cross is a little tricky because it's not an odd size (2x2 rather than 3x3 or 5x5). I've done it using numpy+scipy using a padded 3x3 convolution mask.

import sys
import numpy as np
from scipy import ndimage
import Image

roberts_cross_v = np.array( [[ 0, 0, 0 ],
                             [ 0, 1, 0 ],
                             [ 0, 0,-1 ]] )

roberts_cross_h = np.array( [[ 0, 0, 0 ],
                             [ 0, 0, 1 ],
                             [ 0,-1, 0 ]] )
def load_image( infilename ) :
    img = Image.open( infilename )
    img.load()
    # note signed integer
    return np.asarray( img, dtype="int32" )

def save_image( data, outfilename ) :
    img = Image.fromarray( np.asarray( np.clip(data,0,255), dtype="uint8"), "L" )
    img.save( outfilename )

def roberts_cross( infilename, outfilename ) :
    image = load_image( infilename )

    vertical = ndimage.convolve( image, roberts_cross_v )
    horizontal = ndimage.convolve( image, roberts_cross_h )

    output_image = np.sqrt( np.square(horizontal) + np.square(vertical))

    save_image( output_image, outfilename )

infilename = sys.argv[1]
outfilename = sys.argv[2]
roberts_cross( infilename, outfilename )

From the Wikipedia entry on Robert's Cross. http://en.wikipedia.org/wiki/Roberts_Cross

My script's output.


Why are sidebands generated in AM and FM? When the signal is modulated onto the carrier in the electromagnetic spectrum,
that signal occupies the small portion of the spectrum surrounding the carrier frequency. It also cause sidebands to be generated at frequencies above and below the carrier frequency.

But how and why are those sidebands generated in AM and FM and why are there so many sidebands generated in FM while just two are generated in AM ? Please provide a practical example, as I already know how they are generated mathematically.

What I know is, in the time domain, when the original signal is put into the carrier signal, it is actually multiplied with the carrier signal which means that in frequency-domain the original signal is convolved with the carrier signal. Those two Sidebands in AM are actually the Fourier transform of the carrier signal.

Is this correct?

Carrying information requires bandwidth.

For a given S/N ratio, modulating a signal to carry more information will thus expand its bandwidth.  Call the addition bandwidth "side bands".  If you don't add side bands to a fixed frequency carrier, you can't expand its bandwidth, and thus you can't transmit any information (other than the presence of a constant carrier).

For AM, AM is not PM (phase modulation).  Any additional bandwidth (as required to carry information in the modulating signal) on one side of the carrier will usually have a different phase (change of phase with respect to time from any reference point) from the carrier.  To neutralize this phase difference, AM modulation has to add some additional matching bandwidth on the opposite side of the carrier to carry a signal that will exactly cancel any phase shift of the spectrum the first side, so that AM doesn't become PM.

With FM, modulating a carrier changes the signal frequency to new frequencies.  You can also call those additional new frequencies so generated "side bands".


Precalculate the Normalization of a filter size If we have an Image I and a Gabor filter created at specific size and orientation, assume that is 11x11.

I saw some persons who normalize the size of the filter, but i didn't understand what was meant by "precalculate the normalization of the filter size". And within the normalization, they square the image, for example : Isquare = I^2. So why the image is squared ?

Each pixel of the filter has a magnitude (intensity). The square function equivalently calculate the power density of the pixel. The normalization is to scale the filter power to 1. That's where I^2 at the denominator comes from.


Continuous Time Signal and Discrete Time Signal - Connection Between Periodicity and Discretness I know that all periodic continuous time signal have discrete spectral representations, but are all discrete spectral representations periodic in continuous time?

Also, can all periodic signals be represented by a fourier series?

This is a 2 parts answer:

Given a discrete representation of the Spectral Data does necessarily represent a periodic function in the time domain?
Well, it is by definition. Look, the DFT is a fancy name to a scaled Discrete Fourier Series, namely combination of harmonic signals. By the nature they are all an integer multiplication of a base frequency. Now think about a linear combination of finite harmonic signals which all are multiplication of known frequency -> The are all periodic in the lowest frequency. And so is the signal.
Can all periodic signals be represented by a Fourier Series?
Now you are asking about convergence and there are few requirements about signals to be represented by a Fourier Series which deals with the ability to integrate them over the period. Another issue, when you say Fourier series, remember an infinite number of signals are considered.


Python butter filter: choosing between analog and digital filter types I have a system that performs wireless sampling of data about every 7.5ms (133Hz). Due to it being wireless, I get occasional data drop out. I want to construct a LP butter filter with cut-off frequency of 10Hz using Python's scipy butter method and then downsample everything to a lower frequency.

One of the options that it asks for is whether I want to use analog or digital filter. Isn't running this filter offline in python automatically assume that it's digital? Or does the fact that my samples don't necessarily come in at regular intervals mean that I should stick with analog filter?

I understand how a physical analog filter differs from a digital one, but does this difference apply in the same way to a python-based butter filter?

Isn't running this filter offline in python automatically assume that it's digital?

butter() doesn't filter your signal, it just designs the filter.  It can design an analog filter or a digital filter.

lfilter() is what actually filters your signal, using the filter you designed, and as you can see it is only digital filtering.  It doesn't make sense to filter a digital signal with an analog filter.

If your data is regularly sampled and you want to process it in a computer, then you need a digital filter.  However it sounds like it isn't:

  about every 7.5ms
  ...
  my samples don't necessarily come in at regular intervals

Do you have a timestamp of when each sample was taken? If you do, then do interpolation first:

from scipy import interpolate
f = interpolate.interp1d(timestamps, measurements)
new_timestamps = np.linspace(min(timestamps), max(timestamps), len(timestamps)*3)
new_measurements = f(new_timestamps)

and then digitally filter the interpolated signal. (I'm just picking 3x oversampling arbitrarily.)

Probably your data is not bandlimited? So you need to plot the interpolation and decide which type of interpolation is the most realistic fit for the underlying data.

See Wikipedia: Nonuniform sampling and What is an algorithm to re-sample from a variable rate to a fixed rate?


Camera calibration and extrinsic parameters for perspective transformation I'm working in a photogrammetry camera and laser application where I need to measure a displacement in a specific plane, this is a quick diagram of the setup:

In order to get proper measurements I'm trying to do a projection from the 3d points to the 2d camera plane  proceeding like this:

First I'm doing the Camera calibration using opencv and a chessboard,
I'm taking a few chessboard shots in different angles and applying the
function initUndistortRectifyMap, I´ll have the distortion coefficients,
intrinsic and extrinsic parameters.
With the intrinsic parameters and the coefficients I'll apply lens correction.
The extrinsic parameters of the previous step are for each chessboard
position and I assume that they won't be valid for my particular
application, so I'm placing a chessboard in the above marked area
as "Measurement Plane" and running again the initUndistortRectifyMap function
which this time it's giving me, I think, correct extrinsic parameters
of the ROI.
Once I've got the "correct" intrinsic and extrinsic parameters I´m
applying the camera matrix to calculate the 2d points (u,v) projection:

My questions are:

Is it the correct way to proceed for this application?

The camera system only reads X,Y pixel points each time, but to perform a 2d projection seems like I'll need the Z, how could I have this last coordinate of the 3d point?

For your application, here are the steps I would follow:

Calibrate the camera and obtain the intrinsics - You have obviously done that and corrected for the distortion already. initUndistortRectifyMap only lets you to precompute the projection. Otherwise, you might as well use projectPoints each time. This is a detail and I would not care about that initially. But it's a good practice anyways.)
The main idea is to find the laser plane (with extrinsic calibration to the camera). To find the extrinsics:

a. Place a checkerboard in front of the camera and leave enough space on the board for the laser line to be projected.

b. Image the laser line projections together with the checkerboard.

c. Now we can compute two planes with respect to the camera: The plane joining the laser line with the camera center, and the plane that checkerboard lies in (we now the checkerboard).

d. Intersect the two planes to find the line in 3D space.

e. Collect all such lines from multiple images and fit a plane to all of them. This is the equation of the laser plane.
To make measurements in 3D, simply shoot lines through the observed laser points in 2D image and intersect the rays with the laser plane. Each intersection will give a single point in 3D space.

This way we also do not explicitly need $(R,t)$. Our extrinsics are rather the equation of the laser plane.

In case we use a point-laser, the equation also remains the same, as instead of lines, we can collect sets of points and still do the plane fit.

By the way, here is a toolbox and a recent further read.


How do I compare chords played on a guitar with a known sample? I want to create a project that will recognized the chords that i played in guitar.
I recorded a raw sound to be compare to the one that the user will input.

Example

I have a recorded sound chordA then I will ask the user to play a sound in the guitar and check if the signal is similar to my sound chordA..

Getting the audio signal is easy but comparing them is so hard. What is the method that will suit my application?

This is my code:

recObj = audiorecorder;
disp('Start speaking.')
recordblocking(recObj, 4);
disp('End of Recording.');
myRecording = getaudiodata(recObj);
wavwrite(myRecording,'chordread.wav');
[x, Fs] = wavread('chordread.wav');
[x1, Fs1] = wavread('chordA.wav');
Y = fft(x);
Z = fft(x1);

Now that I have Y and Z how can I compare those two?

You probably need more than one algorithm for an accurate detection. Work on frequency domain (mostly via DFT and spectrogram) is the most often used initial transform. After that sophisticated probabilistic models originally used in speech recognition are applied, such as hidden Markov models, dynamic Bayesian networks, and conditional random fields.

The review article Automatic Chord Transcription Using Computational Models of Musical Context gives you a thorough overview on the relevant technique, although it focuses on a more complicated chords recognition with more than one instrument (key, bass note, beat are all considered as well).

There are also some open source chord detection algorithms:
Chordino, Chordata, and
LabROSA Chord Recognition. Note that in LabROSA, a supervised learning method is applied to "teach" the model recognize the chords with higher accuracy.

This is a very interesting topic, and I look forward to see your update in your guitar chord detection result.


survey/book on adaptive/self-tuning of PI controllers? Can someone point me towards a good book or survey article on adaptive/self-tuning of PI controllers?

All I can find online are articles which are either vague or heavy on theory.

This is definitively not a satisfying answer, but what you call "adaptive" in combination with a pi-controller is commonly called gain-scheduling or lpv-gain-scheduling (lpv = linear parameter varying). Maybe that helps you find a book.


Question about harmonic calculations and other functionals for image analysis I am implementing a method known as trace transform for image analysis. The algorithm extracts a lot of features of an image (in my case features pertained to texture) using a set of transforms on the values of the pixels. Here is a paper describing the algorithm: (PDF link).

The algorithm finally constructs a "triple feature" which is a number that characterizes the image and was calculated by applying a first set of transforms on the pixel values extracted by all the tracing lines going through the image at all angles. On the paper you can see this description starting on page 2, Figure 1 shows an abstraction of such a tracing line which is defined by its angle (phi) and the distance (d) of the line from the center of the image. So we can describe each tracing line by the pair (phi,d).

Now on the next page on the paper (page 3) we have Table 1 which is the set of functionals that is applied to the pixel values extracted by the tracing lines. Each of these functionals applied to all the tracing lines generates another representation of the image. Then another set of functionals (those that appear on Table 2 on page 4) are applied to the new representation of the image along the columns which is basically along parameter d generating, this way, another new representation of the image.

Finally another set of functionals is applied to this last representation which is an array of values over parameter phi which was the angle parameter (that means we have 360 values describing the image, one for each angle of the tracing lines). You can see these functionals in the paper on Table 3 page 4.

Now some of this last set of functionals that have to be applied to the last representation of the image over parameter phi need to calculate the harmonics of this set of 360 values. First I need to calculate the first harmonic and then the amplitude and phase of the first through fourth harmonic but I don't know how to do this.

You can read the first couple of pages of the paper since it's probably better explained there. And figure 1 makes it clear what I mean by tracing lines and the representation of each line by the (phi,d) pair.

Other functionals I have trouble with since I don't know if I'm understanding them correctly are:

$$
 \text{c so that:}\sum_{i=1}^cx_i=\sum_c^Nx_i.
$$

where $x_i$ is one of the values and $N$ is the total number of values.

So my doubt is with what the result should be. The source (the paper listed above) has nothing else to contextualize this functional that's all it says so I'm interpreting that the result is $c$ rather than the sum of $x_i$ until $c$ right?

So let's say I have this set of numbers {2 2 3 1} where each number is one value $x_i$

So the above functional will tell me $c=2$ but the result of the sum is 4. So which should I take as the result in this case given the above functional?

The other functionals I'm having trouble with are related to the calculation of harmonics based on the dataset extracted from the image:

$$
\text{i so that } x_i=\min_{i=1 \ldots N}x_i\text{ without the first Harmonic}
$$
$$
\text{i so that } x_i=\max_{i=1 \ldots N} x_i\text{ without the first Harmonic}
$$
$$
\text{Amplitud of the first Harmonic}
$$
$$
\text{Phase of the first Harmonic}
$$

all the way to

$$
\text{Amplitud of the fourth Harmonic}
$$
$$
\text{Phase of the fourth Harmonic}
$$

I really don't know how to calculate a Harmonic or the amplitude or phase of it for the values. If someone more knowledgeable than me on image analysis can look into the paper linked above and provide me some tips I'll be greatly appreciated.

Perhaps looking at this paper 1 would be appropriate for a literature survey.

The paper describes in detail the rules of constructing functionals for a Trace transform. The functionals described in your paper are merely examples of the functions which follow the rules presented in paper 1.


Understanding this Gabor Filter equation I am working on a project where I need to implement road edge detection. This research paper suits my requirement in the best way, so I am studying it: http://imagine.enpc.fr/publications/papers/TIP10a.pdf

But I cannot understand the Gabor Filter equation they have used and also the set of parameter values they have used. Almost everywhere the Gabor Filter equation is different then theirs. If someone can explain, it would be great.

I will describe the rationale behind the "non-standard" terms, with respect to more classical expressions.

Essentialy, the $\omega e^{-\omega^2 }$ term indicates a derivative of a Gaussian (with 0 average), similar to an antisymmetric wavelet, instead of the modulated Gaussian window in standard Gabor filters: a wavelet filter bank, versus a modulated filter banks. This wavelet "smoothed derivative" aspect might help detect edges.

Then, the $4a^2+b^2$ induces an elliptic shape constraining anisotropy in the 2D wavelet, with a $2$-fold factor ratio.

Last, the $e^{-c^2/2}$ might serve an admissibility condition for this Morlet-Gabor like wavelet, to ensure easier inversion/reconstruction.

It is truly multiscale, as they specify they have 5 dyadic scales.


Negative result from filtering a positive signal with a band-pass filter I have used scipy.signal.remez to calculate the coefficients for a band-pass filter and when I use it to filter a sinusoidal signal that goes between 0 and some positive number (e.g. $2^{16}$) I get that the filtered signal is negative ( it is attenuated, as I expected ) but it has a negative offset, is this expected? I need my outputs to be between 0 and $2^{16}$ (since I am implementing this filter in an FPGA and send the filtered signal to a DAC which accepts positive numbers between 0 and $2^{16}$ Should I just deal with this by adding a positive offset to re-centre the attenuated signal?

Code used to generate the FIR filter:

fs = 2500000         # Sample rate, Hz
band = [95000, 105000]  # Desired pass band, Hz
trans_width = 5000    # Width of transition from pass band to stop band, Hz
numtaps = 200       # Size of the FIR filter.

edges = [0, band[0] - trans_width,
     band[0], band[1],
     band[1] + trans_width, 0.5*fs]
taps = scipy.signal.remez(numtaps, edges, [0, 1, 0], Hz=fs)

The taps are real and look like this:

Here is an example of the result of using this FIR filter to filter a 200KHz sine wave in python:

What you see is what one would expect. As pointed out in Marcus Müller's answer, your band pass filter has a relatively poor stop band attenuation, and since DC is in the filter's stop band, DC is not sufficiently attenuated. You can easily predict what is going to happen:

compute the filter's DC gain, which is just the sum over all filter coefficients: $H(0)=\sum_nh[n]$ This will be a not so small negative number.
Take the DC offset of your input signal ($2^{15}$), and multiply it with $H(0)$. This must be the DC offset of your output signal.

So the DC offset of the output signal can be easily predicted, even without computing the output!

What can you do to get rid the DC offset of the output signal? Since you can perfectly predict it, you can just subtract the DC offset $=2^{15}\cdot H(0)$ from your output signal. Another fun solution is to increase the filter length by $1$ sample, and set that new last sample to the value $-H(0)$ (as computed above). This will make sure that $H(0)$ of the new filter is exactly $0$, and it will only insignificantly change the overall frequency response of the filter.

For the ones who wonder if it is normal that the first and last sample of the impulse response are much larger than all other samples: yes, this is normal, especially for narrow band pass filters, and it is a direct consequence of the equi-ripple optimality criterion implied by the Remez algorithm. These two impulses (first and last sample) are echos generating the sinusoidal (i.e., equi-ripple) stop band behavior in the frequency domain. The ripple in the impulse response generates the (impulse-like) narrow pass band in the frequency domain. So everything can be explained by the basic Fourier relationship "sinusoid $\Longleftrightarrow$ impulse".


What exactly is captured in a Sample of the Nyquist-Theorem variety? I'm new to DSP and was thinking of sampling in the traditional music creation way, in that you capture a sound and play it back, warp it, etc., but when I thought about the Nyquist theorem, I realized that can't be what's happening when sound is being captured.

If we treated a sound clip that is 1 second long as its own waveform, then it's 1Hz and according to Nyquist you'd have to capture it at least at 2Hz. That's only 2 snapshots for 1 second of sound.

What are those 2 snapshots capturing? Let's say it's 16 bits of depth. What can that combination of 1s and 0s capture exactly? I've caught that the bit depth determines how many dBs you have to work with (more range), but I'm a bit lost on the rest. Is it capturing a unique waveform tone? Is it specifically a Sine wave tone? Is it multiple tones? I'm having trouble picturing what's being sampled and how it's being accurately played back.

It might be better to think of Nyquist sampling as shooting fish in a bucket.

The bucket limits what is in the bucket by being only so large. This is like a band limit. We can only have so much variation.

Two samples is the very least number of samples we need to specify that 1 sec of 1 Hz sine wave, although intuition tells us we are on the edge.  Three samples for a second and half is a safer proposition.   How about 120 samples for 60 seconds of 1 Hz signal. If you look at the reconstruction formula for Nyquist sampling, it says you need a need an infinite number of samples to get that 1  Hz sine back.  The summation limits are plus and minus infinity.

In practice, we usually sample a bit faster than what Nyquist specifies.

There are circumstances where one can sample less than what Nyquist tells us and there is a lot of activity in an area known as sparse sampling.

To be honest, this was a terrible analogy. The mathematics is really more straightforward. Things like mathematical limits are hard to analogize.

It also turns out that a 1Hz sine wave is also a mathematical ideal.


Limitations of Arithmetic coding? I have been working on an assignment for university in which I have to compress an 80K message consisting of only two symbols, 'a' and 'b', having probability values of 0.6 and 0.4.

I have successfully written my code in Matlab and it is working efficiently for 10-20 symbols. However when I pass my 80k symbol message through my code it stops updating at the 1057th value due to the Range parameter, which achieves a value of 10^-324 (the theoretical limit of type double float).

My question is, is this the limitation of arithmetic coding when the range becomes so small that we cannot work with it or am doing something wrong?

I am working on arithmetic encoder, so no answers related to range decoding will work.

You should not be using floats or doubles to store your message. You need to write your compressed message in binary.


HPS Algorithm for detecting the fundamental frequency of a guitar string When you plucked on a string, the string is vibrate and except the fundamental frequency,
there are an integer multiple of the fundamental frequency, which are the harmonics.

I know that the HPS algorithm is a tool to find the fundamental frequency, but
I don't understand the downsample part. How can I know how many times I need to downsample
the signal?

I read somewhere that the fundamental frequency can be given by measuring the harmonic frequencies and their GCD (greatest common division)

another source says that there's a need to use the Hanning method.
can someone please explain this algorithm?

OK Lets go!

I consider HPS a very simple algorithm, this image representation show exactly how HPS work !

Yeah is highly recommended apply a hann window for every framed data !

What do you need to do is apply one window function over your framed signal, then apply FFT, you need just the first half absolute values from FFT, its give you the Magnitude, now downsample it for N time, Multiply all downsampled and at  the end find the peak.

  How can I know how many times I need to downsample the signal?

Will depend, the more time you downsample you signal more accurate HPS can be, this has a cost that limits the frequency analysis, for example, imagine that you have a signal sampled at 44100hz, you use FFTsize = 4096 and you decide to make 5 time downsample, its tell me that you can not pick up frequencies above 4403Hz!

floor(4096 / 5) =  819
floor(819 /2) = 409
floor((44100/4096) * 409) = 4403


Magnitude of function in $z$ domain I am newbie to $\mathcal Z$-transform, I searched to find the magnitude of a function in $z$-domain, but I couldn't find anything, for example when we have
$$
H(z) = \frac{z-3}{z-0.5}
$$
How do you compute $\lvert H(z)\rvert^2$? I know that $\lvert H(z)\rvert^2 = H(z)\cdot H\left(z^{-1}\right)$, is it right? when I calculate it for my example:
$$
H(z)\cdot H\left(z^{-1}\right) = \frac{z-3}{z-0.5}\cdot\frac{z^{-1}-3}{z^{-1}-0.5} = \frac{10-3z-3z^{-1}}{1.25-0.5z-0.5z^{-1}}
$$
but it is not a numerical value like Fourier domain.

The $z$-transform can be evaluated at any point on the complex plane that is also in the ROC of the $z$-transform. To find the magnitude of $H(z)$, you can find the magnitude of numerator and denumerator separately, and then divide the results.

Let's say $H(z)=\frac{A(z)}{B(z)}$
So

$$|H(z)|=\frac{|A(z)|}{|B(z)|}$$
and let's assume $A(z)=z-a$ and $B(z)=z-b$ .
$$|A(z)|=|z-a|$$
Let $z=re^{j\phi}$. So
$$\begin{align}
|A(z)|&=|re^{j\phi}-a|\\&=|(r\cos(\phi)-a)+j r\sin(\phi)|\\
&=\sqrt{(r\cos(\phi)-a)^2+r^2\sin^2(\phi)}\\
&=\sqrt{r^2(\sin^2(\phi)+\cos^2(\phi))+a^2-2ar\cos(\phi)}\\
&=\sqrt{r^2+a^2-2ar\cos(\phi)}
\end{align}$$
Similarly $$|B(z)|=\sqrt{r^2+b^2-2br\cos(\phi)}$$
and

$$|H(z)|=\frac{\sqrt{r^2+a^2-2ar\cos(\phi)}}{\sqrt{r^2+b^2-2br\cos(\phi)}}$$


Understanding I and Q modulation As a newbie I am trying to understand I and Q modulation and demodulation. So from what I read the I and Q components of a signal separately mix with RF signals that are 90 degrees out of phase with each other and then the summation is transmitted as an RF signal.

So let's say I am listening to a voice of a reporter via FM radio. The reporter's sound wave gets converted to electricity with a microphone and then presumably gets converted to EM wave with an antenna. I am kind of lost where modulation happens and how would his sound voice gets divided to I and Q before modulation. Is it done mathematically and digitally?

Analog QAM as an example:

https://en.wikipedia.org/wiki/Quadrature_amplitude_modulation#Analog_QAM

same idea for other IQ modulations schemes.


SNR in analog and digital domain In the digital domain a higher sampling will have a better SNR.

Similarly does a narrow filter give a better SNR in analog domain?
Once analog signal is converted into digital signal can SNR be increased? By using a digital filter?

In the digital domain a higher sampling will have a better SNR

No, that's not the case.

In fact, if you oversample a bandlimited signal, you don't get overproportionally more signal power than you get more noise power, but as long as you don't use perfect filters which don't exist, you do see more noise, introduced at every step of the analog-to-digital conversion.

You can get an improved SNR from an oversampled signal. But it's not due to higher sampling, but due to processing gain of what you you do with the sampled signal.

  Similarly does a narrow filter give a better SNR in analog domain?

If you assume your noise has a wider bandwidth than your signal, then quite directly follows that the less bandwidth that contains only noise you let through, the better your SNR gets.

Don't forget that SNR is usually just a "helper figure", not inherently something you care about without looking at what has happened to your signal. For example, inevitably, very narrow analog filters have very nonlinear phase, which leads to very strong distortion of the signal's phase. At the end you might have a very good SNR, which only compares powers of signal to noise, but your signal might still be very, very broken.

  Once analog signal is converted into digital signal can SNR be increased? By using a digital filter?

That is one way of achieving processing gain as outlined above, yes. So, no, your higher-sampled signal doesn't have better SNR, your digitally postprocessed signal has.


2D convolution of image with filter as successive 1D convolutions I want to prove (or more precisely experiment with) the idea that a 2D convoltion as produced by the Matlab conv2() function between an image I (2D matrix) and a kernel (smaller 2D matrix) can be implemented as some 1D conv i.e. the Matlab conv() function and NOT conv2(). Of course possibly some reshapes and matrix multiply might be needed but no conv2().

And to make it clear, I am NOT refering to that kind if thing:

s1=[1,0,-1]'
s2=[1 2 1]
diff=conv2(x,y)-conv2(conv2(x,s1),s2)

diff is = 0 everywhere

Rather, I want to do something like

conv(conv(x(:), filter1)filter2) ...

If a 2D $K_2$ filter kernel is of rank $0$ or $1$, it can be written as a separable product of $2$ 1D kernels   $K_1^r$ and $K_1^c$ on rows and columns. As such, it can implemented by 1D convolutions, as long as one properly reshape the 2D matrices into 1D ones, and take care about "out-of-range" values, to avoid wrap-around. For instance, you can pad in every direction by the size of the filter, and make sure the convolution does not add unwanted information.

Assuming that you know that you have a separable 2D filter, the following code does the job. A one-liner would be:

xRowFull = reshape(conv(reshape(reshape( conv(x(:),s1,'same'),nRow,nCol)',nRow*nCol,1),s2,'same'),nRow,nCol)';

And the code is:

% https://dsp.stackexchange.com/questions/62115/2d-convolution-of-image-with-filter-as-successive-1d-convolutions
%% Initialization
clear all
nRow = 16;
nCol = 16;
HalfSizeCentralImageKernel = 1;
x = zeros(nRow,nCol);
x(nRow/2-HalfSizeCentralImageKernel:nRow/2+HalfSizeCentralImageKernel,nCol/2-HalfSizeCentralImageKernel:nCol/2+HalfSizeCentralImageKernel)=rand(2*HalfSizeCentralImageKernel+1);

%% Original 2D version
s1=[1,0,-1]';
s2=[1 2 1];
y = s1*s2;

%% Step by step 2x1D version
xRowFlat1 = x(:);
xRowFlat1FiltCol = conv(xRowFlat1,s1,'same');
xRowFlat2 = (reshape(xRowFlat1FiltCol,nRow,nCol))';
xRowFlat2 = xRowFlat2(:);
xRowFlat2FiltRowFlat = conv(xRowFlat2,s2,'same');
xRowFlatFilt2Row = reshape(xRowFlat2FiltRowFlat,nRow,nCol)';

%% Compact vectorized 1D version
xRowFull = reshape(conv(reshape(reshape( conv(x(:),s1,'same'),nRow,nCol)',nRow*nCol,1),s2,'same'),nRow,nCol)';

%% Display
figure(1);
imagesc(x);

figure(2);
subplot(1,3,1)
imagesc([conv2(x,y,'same')]); xlabel('Original')
subplot(1,3,2)
imagesc(xRowFlatFilt2Row); xlabel('Separable, step by step')
subplot(1,3,3)
imagesc(xRowFull); xlabel('Separable, one-liner')

diff1=conv2(x,y,'same')-conv2(conv2(x,s1,'same'),s2,'same');
disp(['Max error 1: ',num2str(max(abs(diff1(:))))]);

diff2=conv2(x,y,'same')-xRowFlatFilt2Row;
disp(['Max error 2: ',num2str(max(abs(diff2(:))))]);

[First answer]

Here is a crude  Matlab code. Can you test it, and if OK, i'll send a one-liner (if I can).

nRow = 8;
nCol = 8;
HalfSizeCentralKernel = 1;
x = zeros(nRow,nCol);
x(nRow/2-HalfSizeCentralKernel:nRow/2+HalfSizeCentralKernel,nCol/2-HalfSizeCentralKernel:nCol/2+HalfSizeCentralKernel)=rand(2*HalfSizeCentralKernel+1);
figure(1);
imagesc(x);

% 2D version
s1=[1,0,-1]';
s2=[1 2 1];
y = s1*s2;
diff1=conv2(x,y,'same')-conv2(conv2(x,s1,'same'),s2,'same');
disp(['Max error 1: ',num2str(max(abs(diff1(:))))]);

% 1D version
xRowFlat1 = x(:);
xRowFlat1FiltCol = conv(xRowFlat1,s1,'same');
xRowFlat2 = (reshape(xRowFlat1FiltCol,nRow,nCol))';
xRowFlat2 = xRowFlat2(:);
xRowFlat2FiltRow = conv(xRowFlat2,s2,'same');
xRowFlatFilt2Row = reshape(xRowFlat2FiltRow,nRow,nCol)';

figure(2);
subplot(1,2,1)
imagesc([conv2(x,y,'same')])
subplot(1,2,2)
imagesc(xRowFlatFilt2Row)


Audio processing. Is it possible to directly access the decoded audio data going into the analog input of a computer When a computer transcodes an audio file from one file format to another, does the computer first decode it into the raw digital stream (exactly what is fed into the DAC for the audio output) , or does it do something different without having to completely decode it?

Take Spotify files for example, I believe that they are not in a standard MP3 format that can be played on any device. However, shouldn't it be relatively easy to convert these files into standard MP3 files, by first decoding them and then re-encoding to MP3 format?

Or would it require getting the analogue audio, converting it back to digital and then encoding again. My main query is, do we generally have access to this raw decoded digital audio data?

A program to convert an .mp3 audio file-format into , say, an .ra (real audio) audio file-format needs fully to decode the mp3 file into raw waveform audio and then re-encode it into its new format.

This raw audio waveform data can be contained within 32/64-bit floating point or some integer formats though. But when it's sent to audio DAC, it should be in an integer format; this depends on the PC audio standard and as far as I know, Windows (and like) systems use 8/16 bit integer formats at this stage (you have to check for their developer knowledge base for their latest implementations)

So there can be a very slight difference between the raw data at the input buffer of DAC and raw data that reside in the system RAM, but they are both raw waveform data.

Nevertheless, unless the the two audio formats have some common design layering, then it's not possible to transcode from one format to the other without full decoding of the original.

Finally, analog audio is practically never required as it would introduce more errors and complexities than do anything useful.


Filtering of streaming audio Given something like 103 data points ($N=103$), a DFT will give back 103 frequency values. Then to do something like filtering the high frequencies involves setting the high frequency values from the DFT to zero, and doing the inverse DFT to get back 103 data points that represent the original signal without the high frequencies.

This makes sense to me given all 103 data points all at once. But what about streaming an audio WAV file that is rather large (say for example $N=10^5$). If one wanted to filter the high frequencies then the approach i just described on the entire data of $10^5$ points logically makes sense. But that is not reasonable when streaming a WAV file for playback. What is done to high frequency filter a streaming playback of an audio file?

FFT -> zeroing coefficients -> IFFT is not the correct way of doing filtering - the actual filter realized by doing so has poor characteristics.

The correct way of filtering signals is to compute the coefficients of a digital filter, a process known as filter design and for which a large body of software tools/documentation is available, and apply it to your input sequence. In short, this consists in evaluating for each sample a linear combination of the past input samples, and the past output samples. Depending on the requirements of your filter in terms of stop-band rejection / ripple, only a few coefficients could be necessary, making it way more efficient than FFT. Since the only information required to compute an output sample is the few past input/output samples, there is no problem to apply it to streaming audio.

You will need to use FFT only if you decide to go with FIR filters, and if your filter requirements cause them to have an insanely large number of coefficients. In this particular case, it will be efficient to apply the filter to successive blocks of your input data through FFT and overlap-add.


Non-cyclic smoothing of a 2D image I was wondering if there was a simple approach to smooth images, e.g using a Gaussian Kernel without introducing boundary effects which are also present when using Discrete Fourier Transforms?
For both approaches, the missing pixels at the left-hand side will be taken from the right-hand side and vice versa.

Right now I have a 5x5 Gaussian Kernel which I move iteratively over my MxN image. Up to now, in case I reach the last pixel of the image, the operator automatically uses  the first pixel and so on. Since my picture is not periodic I get artifacts even if I pad my picture on all sides with the intensity values of the corresponding adjacent column or row.

Do you gus have any suggestions?

Thanks!

Pad image and smooth:

Make image of value 1 same size as original. Pad and smooth:

Divide original image by '1' image:

And crop:

MATLAB demo code:

I = double(imread('lena.tif'));
g = fspecial('gaussian',60,10);
Ipad = padarray(I,[10,10]);
imagesc(Ipad); colormap gray(256); pause;
Ipadg = imfilter(Ipad,g);
imagesc(Ipadg); colormap gray(256); pause;
A = ones(size(I));
Apad = padarray(A,[10,10]);
imagesc(Apad); colormap gray(256); pause;
Apadg = imfilter(Apad,g);
imagesc(Apadg); colormap gray(256); pause;
IA = Ipadg./Apadg;
imagesc(IA); colormap gray(256); pause;
I = IA(11:end-10,11:end-10);
imagesc(I); colormap gray(256); pause;

*Note: * MATLAB imfilter already does zero padding. The pad array step is not needed, it is just for demonstration. Actual code would be:

I = double(imread('lena.tif'));
g = fspecial('gaussian',60,10);
Ig = imfilter(I,g);
A = ones(size(I));
Ag = imfilter(A,g);;
Ismooth = Ig./Ag;


Fourier transform of $\sum_{n=-\infty }^\infty e^{-|t-2n|}$ I need to find Fourier transform of

$$\sum_{n=-\infty }^\infty  e^{-\lvert t-2n\rvert }$$

I found the solution online. But I think the choice of limits in equation (1) of the solution is wrong. It should be probably from $-\infty$ to $2n$ for the first integral and from $2n$ to $\infty$ for the second integral. But I am confused about what the limits would be as the both integration and summation is involved. and I am not sure whose limits to change (probably limits of integration must be changed only and that of summation must be kept as it is)

Given Solution:

Please tell me what the limits will be and the explanation. Also if possible please solve the problem.

Original answer:

$$ \frac{(1-e^{-2(1+jw)})}{(1-e^{-2})(1+jw)}-\frac{(e^{-2})(1-e^{-2(1+jw)})}{(1-e^{-2})(1-jw)}$$

First of all, the signal is periodic, so its spectrum must be discrete. Hence, your final solution ("Original answer") cannot be correct.

I would go the following way:

Calculate $X(f)=\mathcal{F}\left\{\exp(-\lvert t\rvert)\right\}$
Apply the time-shifting property of the Fourier Transform to transform all time-shifts of the function.
See that you can factor out $X(f)$ from the resulting sum and then see that the terms in the sum correspond to a comb function where the diracs are $\frac{1}{2}$ apart from each other.

Edit: Below is some code to plot the function:

Fs = 100  # some sampling frequency for discretizing the signal to plot it
T = 10    # the time range
t = np.arange(-T, T, 1/Fs)  # the time samples

# define the basis function
f = lambda t: np.exp(-abs(t))

# plot the basis function and the periodic extension
plt.plot(t, f(t))
plt.plot(t, sum(f(t-2*n) for n in range(-10, 10)))

Edit: Regarding your comment:

$$
\begin{align}
\mathcal{F}\left\{e^{−\lvert t\rvert}\right\}&=X_1(j\omega)=\frac{2}{1+\omega^2}\\
\mathcal{F}\left\{\sum_{n}e^{-\lvert t-2n\rvert}\right\}&=\sum_nX(j\omega)\exp\left(j2n \omega\right)\\
&=X(j\omega)\sum_n\exp(j2n \omega)
\end{align}
$$

Here's the trick:
$$\sum_n\exp(j2n \omega)=\frac{1}{2}\sum_n\delta\left(\omega-\frac{2\pi}{2}n\right)$$

Have a look at the wikipedia article and also this article for example on why this holds.


Blind Signal separation vs .... what? I've been learning about blind signal separation.

I have noticed that there are a number of algorithms for performing blind signal separation. (eg: Independent component analysis)

Which begs the question, are there any 'non-blind' signal separation algorithms?

I've been searching for a while and I haven't found anything. Nothing is listed on the wiki page and I've never noticed any in open libraries such as open CV.

I'm not even sure what search terms to look for.

The opposite of blind-signal separation is not "non-blind signal separation". The thing is, depending on what you know on the signals you will use different algorithms.

For example, if you know that there is a low frequency signal mixed with a higher frequency one, you can use low/high-pass filtering to get the two signals.

There's no real name for what you're looking for because every case is different. The generic way is more or less filtering the signals with the information you know. But the filtering is always different.


How to design a variable slope highpass filter I want to design some specific slope highpass filter to enhance  speech; just like $1$ dB/octave,  $2$ dB/octave and more. Is there a method to calculate the coefficients by specific slope?

The order of a filter determines its slope. Normally orders are positive integers (1, 2, 3, ...). So the slopes are fixed to multiples of approximately 3db/octave. However, taking the case of a Butterworth filter, someone has applied partial calculus to obtain fractional order Butterworth filters which will effectively allow you to chose the slope. If you do a Google search for fractional order Butterworth filters you will get back some papers on the subject. For instance here: https://arxiv.org/ftp/arxiv/papers/1210/1210.8194.pdf

If you need a different filter type you could apply the same thinking, but the maths may be more complicated for, say, an elliptic filter.


Continuous frequency vs discrete frequency? I have understood idea of discrete time and continuous time but I am feeling difficult to comprehend this idea in regard to frequency.

As for example the DFT output is discrete and the DTFT output is continuous.

How we can realize/visualize difference between continuous frequency and discrete frequency?

Both continuous-time and discrete-time signals generally have a continuous spectrum. Only periodic signal have a discrete spectrum. You might be familiar with the Fourier series, which is just a frequency-discrete representation of a periodic signal. There is also a Fourier series representation of discrete-time periodic signals, which is of course also discrete in frequency.

The discrete Fourier transform (DFT) is basically a Fourier series representation of a finite length discrete-time signal, which is thought of as periodically continued outside its support. The actual spectrum of the finite-length signal is continuous in frequency, and the DFT computes equidistant samples of this frequency-continuous spectrum.

The DFT can be used to compute discrete approximations to the continuous-time as well as the discrete-time Fourier transforms, both of which are continuous in frequency. The DFT is used so frequently, because we have efficient algorithms for its computation. These algorithms are called Fast Fourier Transform (FFT) algorithms.


Anti-aliased boolean masks for integration of CCD images (This question was already asked on Stack Overflow under the tags Matlab/Octave/Signal-processing/Image-processing, but I am reposting it here following the advice of one user that this site may be a more appropriate forum for such questions.)

Is it common/useful to generate an anti-aliased boolean mask for CCD image integration purposes ?

Let's imagine that I want to integrate the light that strikes a given area of a CCD sensor. The boundaries of this area correspond to physical coordinates that have no obligation to coincide exactly with the pixels of my CCD sensor.

In Matlab/Octave/(or even Scientific Python), the common algorithm used for integrating such an area is to define a boolean mask using a logical operation on an array, such as:

mask = R < radius;

(where 'R', and therefore 'mask', are 2D arrays and 'radius' is a float). Such a mask has values of 0 or 1.

I can then integrate the pixels that are comprised within the boundaries of my area by summing the masked image:

integrated_signal = sum(sum(mask.*image));

(where 'image' is the output of the CCD sensor, a 2D array with U16 values for instance).

However, mathematically speaking, nothing would prevent me from defining a continuous-valued mask, that is, a mask whose pixels could take any value between 0 and 1. Physically, this means that I could integrate portions of pixels. I could even compute such a mask using an anti-aliasing algorithm in order to create better approximations of non-square-friendly shapes such as circles (for integrating a signal with revolution symmetry, let's say).

My question to digital imaging geeks is the following: Is it standard practice to define such anti-aliased mask for integration purposes of imaging sensors ? Is it useful (as I suspect), or does it produce physically irrelevant results ?

You're right: mathematically speaking, nothing would prevent you from defining a grayscale mask. Practically, it is not a good approximation of what is going on.

A mask value of 0.5 means that you'd count half the pixel value to your signal. First of all, this will require you to exactly estimate and subtract the camera offset, b/c adding half the offset is going to screw up your result. More importantly, this assumes that half of the signal value on your pixel must have come from somewhere else.

There are three sources for where your other 0.5 of the signal may have come from:

If the source whose signal you're integrating is isolated, all the signal on that edge pixel will have come from your source of interest. So multiplying that by 0.5 is wrong.
If there is a second source that also might contribute to the pixel, the proper ratio for dividing up the signal depends on the relative strengths of each signal, and on the relative positions of the center. While a grayscale mask makes sense here, you have to calculate it iteratively, since it depends on the signal integral.
If there is noise on the sensor, you'd be claiming that the other 0.5 is noise. That's most likely wrong as well - you want to add the full noise to your signal so it averages out properly (that can be a bit of an issue with #2).

This is why I have never encountered grayscale masks except for case #2. In addition to being a bad approximation to reality, grayscale masks are less convenient to handle. With logical masks, you can, for example, conveniently write

integrated_signal = sum(image(mask));


Encoding sequence of unfair coin flips The question: Consider transmitting the results of $1000$ flips of an unfair coin where the probability of heads is given by $p_H$. The information contained in an unfair coin flip can be computed:
$p_H\log_{2}(1/p_H)+(1−p_H)\log_{2}(1/(1−p_H))$
For $pH=0.999$, this entropy evaluates to $.0114.$ Can you think of a way to encode $1000$ unfair coin flips using, on average, just $11.4$ bits? (question from https://web.mit.edu/6.02/www/f2011/handouts/2.pdf)
My wrong answer: I thought that I could encode the location of the bits which turn up tails. Since there are 1000 flips, I could encode every flip using 10 bits ($2^{10}=1024$). taking the average expected length to encode each flip and then multiplying by $1000$ for all the flips gives:
$1000[(0.999)(0)+(0.001)(10)]\\
1000(0.001)(10)\\
10$
But I know that any encoding which averages a smaller length in bits than the entropy must have some ambiguity in the message, so since $10<11.4$, what information is my coding system missing?

The problem is the assumed knowledge that the receiver needs to have. In your coding scheme you assume that the receiver knows that you transmit exactly $1000$ symbols. If the receiver didn't know that, there is no way to distinguish the following two cases:

$2$ tails at certain positions inside one block of $1000$ symbols
$2$ tails at the same positions as for case 1. but in two different blocks

So your coding scheme is incomplete and that's why you end up with an average bit rate less than what we would expect from the entropy of the source.
Also note that your scheme, even though incomplete, will exceed the minimum possible bit rate for larger block lengths. That limit of the block length beyond which the required rate becomes larger than the minimum possible rate can be computed as the smallest integer value of $N$ satisfying
$$p\lceil\log_2(N)\rceil>-p\log_2(p)-(1-p)\log_2(1-p)\tag{1}$$
where $p$ denotes the probability of a tail.
Apart from the above, you would also need to assign a codeword to the case that there is no tail inside a given block. Of course, for this likely case it would be wise to choose a short codeword.


Solving for frequency deviation, $\Delta f $ (FM modulation) I'm sorry if I posted this in the wrong section. I solved the question and got $50\textrm{kHz}$, but the solution says it is $100\textrm{kHz}$ without any proof.

Im asking for someone to tell me which is correct.

If you use the standard definition of frequency deviation as the

  "maximum instantaneous difference between an FM modulated frequency and the nominal carrier frequency"

then for the given example you indeed have $\Delta f=50\;\text{kHz}$, because the nominal carrier frequency is $f_c=10\;\text{MHz}$, the maximum frequency is $f_c+\Delta f=10.05\;\text{MHz}$, and the minimum frequency is $f_c-\Delta f=9.95\;\text{MHz}$.


How to find the period of a noisy signal using MATLAB's $\tt xcorr$? Here is my code:

function [ T ] = FindPeriodicity2(x,Fs)
ac=xcorr(x,x);
[~,locs]=findpeaks(ac);
T=mean(diff(locs)/Fs);
end

and when I pass the signal x with a Gaussian noise the period is 1.564 but as we know it equals $\pi$, so the answer is about half of the real period in this case. I know I should eliminate some of the peaks found in 'ac', but I have no idea how to choose from peaks such that it works for any signals.The plot of 'ac' is attached .
here is the signal x:

Fs=100;
t=0:(1/Fs):30;
n=length(t);
noise=randn(1,n).*0.3;
x=sin(2.*t+pi)+cos(4.*t)+noise;
T=FindPeriodicity2(x,Fs);
disp(T);

any hint would be appreciated.

You are getting that result because you are taking into account other maxima that do not correspond to the points you are actually looking for.

When you call findpeaks() with no other parameter than the input signal, you are taking into account all the points marked with a blue arrow. To determine the period, you just want those inside red circles.

Note that, due to the noise, inside some red circles there are more than just one point considered as a peak. That is undesired as well.

You have to put a restriction to find peaks such that they are at least 160 samples away from each other (in this case, due to your signal and its correlation). You can see that if you replace

[~,locs]=findpeaks(ac);

with

[pks,locs] = findpeaks(ac,'MinPeakDistance',300); %Or any number >160

the period will be the expected one.


Black box signal prediction I have a stateful black box with four real-valued inputs and one real-valued output.  My problem is to predict the output at each moment in time, given the sequence of inputs seen up to that point.  During a learning phase, I can vary the inputs however I want and observe the output.  There's a little noise, of course, and the black box doesn't appear to be completely deterministic.

Specifically, I'm modeling a hard drive, and I want to predict the access time of the latest request given all previous requests.  I want a more black-box approach, though, because of the complexity of explicit models, and because I want this to work for other similar devices such as SSDs.

A couple of people have suggested that signal processing might be appropriate to analyze the sequences of input and output values.

Are there any ideas from signal processing that could help me to predict the output, or to characterize the input?

Generally, for non-linear systems the aren't any tools that are guaranteed to work. You need to know something about the nature of the box. If you can model it with a system with unknown parameters, "learning" by observing input-output relationships can help you estimate those parameters, but I doubt that you can blindly "learn" the system model, especially if it has memory/states. Having said that, a more or less generic method of estimating nonlinear systems up to a polynomial degree is using Volterra kernels with some sort of gradient descent method for recursive error minimization. Methods like LMS and RLS are widely used.

Volterra kernels help you estimate a system of the form

$$y(t) = k_0 + \sum_{n=1}^{\infty} \int_{-\infty}^{\infty}\cdots \int_{-\infty}^{\infty}k_n(t_1,t_2,\ldots,t_n)x(t-t_1)x(t-t_2)\cdots x(t-t_n)dt_1dt_2\cdots dt_n$$

One way to look at the following system is to notice its similarity to convolution with addition of non-linear terms in $x$. Note, however, that if non-linearity of your system cannot be modeled with a polynomial (exponential, logarithmic and many other functions) this will require infinity complexity to estimate your system correctly.

There aren't many papers on the topic available online for free, but you can look at this one and this one to get an idea of what this is all about.


Derive DFT of $x[2n]$ Question is to find DFT, $G_7(k)$ of $g_7(n)$ = $x(2n)$.
For purpose of illustration only, we can assume that $x(n)$ can be represented by envelope shown in corresponding figure.

Shown in Figure Below, is solution to this problem.

But, I am unable to understand this Solution.
If anyone can help me understand this solution, I will be very thankful.

You haven't specified what part exactly of the solution you didn't understand.

The purpose of the posed question is to show the spectral effects of signal downsampling, in this case by a factor of 2.

The original discrete signal is $x[n]$ (whose DTFT is $X(e^{j\omega})$), with a length of N samples, and the downsampled signal (by a factor of 2) is $g_7[n]=x[2n]$, so $g_7[n]$ is made from the even samples of $x[n]$, with a total length of N/2 samples.

The first experssion of the solution is simply the DFT (by definition) of $g_7[n]$. The solution then proceeds to some series manipulation, then uses the DTFT-derived definition of DFT to show the final expression.

The final expression of the solution shows that downsampling the signal by a factor of 2 produces a double-frequency image of the original signal in the frequency domain. This alludes to the importance of anti-aliasing filtering prior to signal downsampling.


Register 2D image - 3D model I have 2D pictures of an airplane of a specific model, taken with a simple camera, and a 3D point cloud model of this airplane.
I would like to project the pictures onto the 3D model (cloud of points) based on control points that I manually locate on both ends. This should allow me to calculate distances on the 2D pictures.
Is there a programming library that does this? I program in Python mostly.

You could easily use OpenCV (https://opencv.org/) to do that. Check out the solvePnP functionality. Here is a demo on head pose estimation.


How to know the filter order Could anyone explain to me how can I know the filter order based on algorithm?

For example:

$y[n] = \frac12 x[n] – x[n–1] + \frac12 x[n–2]$ has order filter 2.
$y[n] =2 x[n] – x[n–1] + y[n–1]$ has order filter 1.
$y[n] =2 x[n] – x[n–1] + x[n–2]+ y[n–1] + y[n–2]$ has order filter 2.

Why?

If the discrete-time filter is described by the linear difference equation

$$y[n]=\sum_{k=1}^{K}a_ky[n-k]+\sum_{m=0}^{M}b_mx[n-m]$$

(with $a_K\neq 0$ and $b_M\neq 0$), the filter order is given by
$\max\{K,M\}$, i.e., by the maximum delay necessary to implement the filter.


Is RMS of a sound wave good measure of loudness? I'm not talking about FFT, just the raw sound wave and the RMS of it.

Thanks

No

The sensitivity of human auditory is very dependent on frequency and without any type of frequency weighting, a raw RMS estimate will be pretty much useless.

There are some standard frequency weighting curves, most popular being the "A" weighting.  See https://en.wikipedia.org/wiki/A-weighting. The sensitivity is also dependent on the absolute level and hence there are multiple weighting curves.

Sp frequency weighted RMS is a good starting points, but they typical ignore any sensitivity to time domain behavior like impulse, transient, masking, etc. There are good models to predict actual perceived loudness but they tend to be rather complicated.


Understanding Ergodicity and Ensemble Averaging Literature says that a stationary signal is ergodic, if its ensemble average = time averages. Should it be the statistics computed by time averaging = statistics computed by ensemble averaging?The way I understood from the book in the link is as follows :

Case1: There are N samples in a time series I compute the mean1,skewness1,kurtosis1 .

Case2: Next, with the N samples, T experiments are conducted and in each of the experiment is calculate the statistics. After, T experiments are done, we calculate the average of these statistics - mean_avg,skewness_avg, kurtosis_avg. This is called as Ensemble averaging or Probabilistic averaging.

If mean1 = mean_avg and so on, the the process/signal is ergodic. Is my understanding correct?

I will try to explain this practically using MATLAB notation.
Yet before that I must say the ergodic property sometime is limited to a level of moment, namely ergodic in the 1st , 2nd, 3rd moment, etc...
If the process is IID it is promised to be ergodic.

Now, assume we have a function myRandomProcess in MATLAB.
It returns a row vector of length n with the random properties of our random process.

Now, assume we build a matrix A which is $ m $ by $ n $ matrix which each row is a different vector created using the above function.

Now, using the mean and var functions in MATLAB we can compute the properties of the process.

If we apply mean a long the columns (Namely, on each row) we will get $ m $ results.
If those are close enough to the $ n $ results yielded by applying the mean function along the rows (On each column by itself) we can assume the process is ergodic in its first moment.

Now, if $ m $ and $ n $ approaches infinity "close enough" can be small as we want.

Most properties of random process related to its ergodicity can be inferred by looking on its auto correlation function.

Update

Simple MATLAB code to convey the idea:

% Ergodic Random Process

% The Prbability dimension is a long rows (Vertical) and the time dimension
% is a long the columns (Horizontal).

numSamples  = 1000;
numTrials   = 1000;

% Ergodic Process - IID Random Samples

mRandomProcess = randn(numTrials, numSamples);

vTimeAverage = mean(mRandomProcess, 2);
vSampleAverage = mean(mRandomProcess, 1);

vTimeVariance = var(mRandomProcess, 1, 2);
vSampleVariance = var(mRandomProcess, 1, 1);

disp(['Ergodic Process Statistics']);
disp(['1st Moment Ergodicity - ', num2str(mean([vTimeAverage(:) - vSampleAverage(:)] .^ 2))]);
disp(['2nd Moment Ergodicity - ', num2str(mean([vTimeVariance(:) - vSampleVariance(:)] .^ 2))]);

% Un Ergodic Process

% The DC Level is drawed once per random process (Realization)
vDcLevel = (10 * rand(numTrials, 1)) - 5;

mRandomProcess = randn(numTrials, numSamples) + repmat(vDcLevel, [1, numSamples]);

vTimeAverage = mean(mRandomProcess, 2);
vSampleAverage = mean(mRandomProcess, 1);

vTimeVariance = var(mRandomProcess, 1, 2);
vSampleVariance = var(mRandomProcess, 1, 1);

disp(['Un Ergodic Process Statistics']);
disp(['1st Moment Ergodicity - ', num2str(mean([vTimeAverage(:) - vSampleAverage(:)] .^ 2))]);
disp(['2nd Moment Ergodicity - ', num2str(mean([vTimeVariance(:) - vSampleVariance(:)] .^ 2))]);

As you'll be able to see, the first process is Ergodic (As any IID process), as it is equivalent to calculate its moment along the time dimension or the sample dimension.

On the other hand, the second one isn't Ergodic.
Since its mean along the rows is zero (Drawn from a uniformly random variable [-5, 5]) yet along its columns its mean is the value of the specific realization of the uniform variable.

I hope this clarifies it.


Image processing using conv-neural net Okay I am starting to understand the idea of neural networks but I still haven't been able to understand the use/benefits/implementation of convolutional neural networks especially for image processing.

What I don't understand is how to detect an object of different size.

Say we have a weights matrix of 50*50 and a face in the training set which is of 50*50 size . Now when you take the weighted sum of the matrix it will return a particular value say "X" . So now I understand the idea is to run the weights matrix over the entire image to get a lock over a region which gives a weighted sum of ~"X" now you have detected a face.

If the above mentioned understanding is true Then how is that an image of size 25*25 going to churn out a value anywhere clone to "X" . So the basic question is how to deal with relative sizes of the image .

If possible please suggest some good tutorials for the same .

First of all, CNNs are typically used with identical sized inputs. So, you resize your input to a constant size, no matter what.

Next, detecting an object and recognizing an object over many samples are different tasks. For object detection, you could be better of with R-CNNs. Its sources are also available and there is a Caffe tutorial for R-CNNs. Of course there are many other works on the same topic. For example, this one.


Algorithm to detect if a spoken phrase matches a previously spoken phrase I would like to create a very simple system that allows a user to train (once or multiple times) a spoken phrase (from 1 word to a whole sentence.)

Then the (same) user would speak the phrase back, and the system would score how close to the original the phrase was matched. If certain words were missed or incorrect, the score would be deducted.

I've seen several articles about audio fingerprinting (ie Shazaam) but these implementations are an overkill. I have a known sound pattern that I'm matching against, not search against a database of audio.

Here's what I've tried (using the R libraries DTW and tuneR):

recorded the phrase "Known as the father of the constitution, he was the first to arrive in Philadelphia for the Constitutional Convention." three times, however, the 3rd time instead of saying "Philadelphia" I said "New York".
took the MFCC of all three WAV files
took the DTW three times: Phili1 vs Phili2, Phili1 vs. NY, Phili2 versus NY
Plotted the alignment, question is...what do I do with this data, what metric should I use You can see in the plot that the first alignment is most linear (same phrase repeated exactly)

Phili1 vs Phili2

Phili1 vs. NY

Phili2 vs. NY

MFCC is an established way of handling speech recognition. You will find plenty of libraries and tools (e.g. yaafe, SoX) that can extract them for you. Implementing them yourself is also quite easy. Yaafe provides plethora of other audio features you might like to play with.

Mind that MFCC are not very robust against noise, so you might like to apply some noise reduction beforehand if required. Some people advocate raising the log-mel-amplitudes to a power 2 or 3 before applying discreet cosine transform.

Some libraries, like scikit learn, next to all machine learning algorithms provide also DTW. On your place I would just try MFCC and see how it goes from there.


How to make a silence? I am trying to make a method in my simple application which is like SetSilence(); Now If I have audio (mono) PCM 16 bit Signed - should I fill all data with 0 values? And for PCM 16 Unsigned with 32768 values?

Yes, this is correct. Other values would cause a DC offset.


Optimally approximating the sign function by functions with compactly supported Fourier transform I'm re-posting a question of mine from math.stackexchange in hopes that folks here might have the right kind of expertise.
I'm looking for a systematic way to approximate the sign function
$$\operatorname{sgn}(x) \triangleq \begin{cases}1&\text{ if }x > 0\\
0&\text{ if }x=0\\
-1&\text{ if }x<0\\\end{cases}$$
by functions $f$ whose (distributional) Fourier transforms $\hat{f}$ are compactly supported.
More precisely, let us fix some $\epsilon$ small. Then I am interested in odd, continuous $f$ such that:

$\text{supp}{\hat{f}}\subseteq[-\pi,\pi];$
there exists some $a_f>0$ such that $|f^2(x)-1|<\epsilon$ whenever $|x|\geq a_f$.

Goal: I am looking for a "best" $f$ in the sense $a_f$ is as small as possible.
I figured that one way to do this would to be to look at convolutions $f=a*sgn$ where $a$ is a function that satisfies $\int_\mathbb{R} a(s)\,ds=1$ and $\text{supp }{\hat{a}}\subseteq[-\pi,\pi]$. For example, we could take

$a_1(s)=\dfrac{\sin(\pi s)}{\pi s}$; or
$a_2(s)=\frac{1}{2}\left(\dfrac{\sin(\frac{\pi}{2}s)}{\frac{\pi}{2}s}\right)^2,$

where the Fourier transforms of $a_1$ and $a_2$ are a rectangular and a triangle pulse respectively.
Doing some numerical experiments, I find that for the choice $\epsilon=0.05$, $a_2$ is "better" than $a_1$ by the above criterion.
Question: Is there a more systematic way to find better (or a best) $f$? References that deal with this type of problem would also be appreciated.

This problem is in a way the dual of the problem of approximating a Hilbert transformer by a filter with finite memory. The frequency response of an ideal Hilbert transformer is
$$H(\omega)=-j\,\textrm{sgn}(\omega)\tag{1}$$
and the corresponding ideal impulse response is
$$h(t)=\frac{1}{\pi t}\tag{2}$$
We can approximate $(1)$ and $(2)$ by multiplying the ideal impulse response by a window of finite length:
$$\tilde{h}(t)=h(t)w(t)\tag{3}$$
where the window function $w(t)$ is a real-valued even and smooth function with finite support.
In your case you could use a window function $W(\omega)$ in the frequency domain to approximate a sign function:
$$f(t)=\textrm{p.v.}\frac{1}{j\pi}\int_{-\pi}^{\pi}\frac{W(\omega)}{\omega}e^{j\omega t}d\omega=\textrm{p.v.}\frac{1}{\pi}\int_{-\pi}^{\pi}\frac{W(\omega)}{\omega}\sin(\omega t)d\omega\tag{4}$$
Note that this is just a generalization of what you already did. You used a rectangle and a triangle for $W(\omega)$. However, there are many more "nicer" window functions that are smooth, and that may result in better approximations of the sign function.


Seamless audio loops via Fourier transforms? I am trying to create short seamless loops of continuous sound.

I recorded myself making an "Aaaaah" sound at 44.1 KHz, and cut this recording to a section containing 4096 samples (giving me slightly under 1/10th of a second of audio).

When I play this in a loop, the sound is very much recognizable but there is of  course an audible ticking noise at the point where it is cut.

I applied an FFT transform to this data. When I apply an inverse FFT to the Fourier tables to recreate the signal, the defect is of course still audible.

Is there anything I can do to the Fourier tables to eliminate the defect without changing the other audible properties of the recording? I'm very much uneducated in signal processing, but intuitively I feel that since the ticking sound repeats once per loop, it might be possible to work in the frequency domain to achieve my goal of cleaning it up to produce a seamless continuous sound.

Note: In case there are Javascript programmers here, I created a CodePen here with the "Aaah" loop and its conversion to and from Fourier tables, to facilitate experimentation in code.

Try using some sort of pitch estimator to determine the pitch of your sample.  Or hum in tune to a note of a known musical pitch (from a pitch pipe for instance).  Then overlap and cross-fade frames, not 4k apart, but a shorter distance apart that is an integer multiple of the period corresponding to that pitch.


Signal processing research idea I am a graduate student in applied math, specializing mostly in harmonic analysis. My undergraduate studies were mostly in pure math but I realized I prefer the applied/computational side. I was hoping to undertake a reasonable project which will lead to a short, but publishable paper in the area of signal processing. I have toyed around with a few ideas but I just cannot come up with something concrete as of yet.

Topics of interest:

Data compression
Distortion/noise
Optimization problems
Image processing

This is very simple. You need an inter-disciplinary domain where all (and more) of your interests are satisfied, while being very publishable, employable and fun. Easy!! Jump into the bleeding edge of Speech Recognition.

There is TONS of current research is Automatic Speech Recognition (ASR). The domain is very multi-disciplinary (applied math, linguistics, computer science, optimization, electrical engineering) with applications to just about every part of life where a speech interface can work (healthcare, mobile assistants, in-car, entertainment/gaming, education, etc)

For example, state of the art systems like Baidu's Deep Speech 2 use convolutional neural networks (CNNs) to process speech spectrograms (DSP!). CNNs have revolutionized image processing since 2012. The CTC objective that is optimized during training is a very interesting loss function for supervised training of sequences.

Finally, have a look at these 2016 Stanford University student projects, notably those on speech recognition. This should give you some ideas of what scope or size of project makes sense.

Conditioning WaveNet on Learned Formant Characterizations for Speech
Audio Enhancement
Music-Speech Discrimination
Semantic Analysis of Political Speeches
Language Recognition and using Speech Recordings
Reinforcement Learning for Feature Selection in Affective Speech Classification

Good luck!!


Different approaches for partial image derivation I know there are different ways for partial derivation of an image, among others: Sobel kernel, LoG, Prewitt and so on.

But the simplest one is the central difference:

$$
  \frac{d}{dx} f(x) \approx \frac{f(x+1) - f(x-1)}{2} \longrightarrow 0.5[1\ 0\ -1]
$$

Which means convolving the image with above matrix.

Assume the image looks like this:

$$
        I(x,y) =
        \left(
        \begin{matrix}
        0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
        0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
        0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 \\
        0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 \\
        0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 \\
        0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 \\
        0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
        0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\
        \end{matrix}
       \right) \in \mathbb R^{8\times8}
$$

Convolving this image with matrix $$G_x = 0.5\ [1\ 0\ -1] \in \mathbb R^{1 \times 3} $$ results in:

$$
\lvert I \ast G \rvert =
        \left(
        \begin{matrix}
 0 & 0 &  0 & 0 & 0  & 0 & 0 & 0 \\
 0 & 0 &  0 & 0 & 0  & 0 & 0 & 0 \\
 0 & 0.5 & 0.5 & 0 & 0 & 0.5 & 0.5 & 0 \\
 0 & 0.5 & 0.5 & 0 & 0 & 0.5 & 0.5 & 0 \\
 0 & 0.5 & 0.5 & 0 & 0 & 0.5 & 0.5 & 0 \\
 0 & 0.5 & 0.5 & 0 & 0 & 0.5 & 0.5 & 0 \\
 0 & 0 &  0 & 0 & 0  & 0 & 0 & 0 \\
 0 & 0 &  0 & 0 & 0  & 0 & 0 & 0 \\
        \end{matrix}
       \right)
$$

The cells containing $0.5$ are the edges of the image in $x$ direction.

Know assume we would extend our filter $G_x$:

$$
G_x = 0.5
\left(
\begin{matrix}
 1 & 0 & -1 \\
 1 & 0 & -1  \\
 1 & 0 & -1 \\
\end{matrix}
\right) \in \mathbb R^{3 \times 3}
$$

Now convolving this filter with our image results to:

$$
\lvert I \ast G \rvert =
        \left(
        \begin{matrix}
 0 & 0 & 0 & 0 & 0  & 0 & 0 & 0 \\
 0 & 1 & 1 & 0 & 0  & 1 & 1 & 0 \\
 0 & 2 & 2 & 0 & 0 & 2 & 2 & 0 \\
 0 & 3 & 3 & 0 & 0 & 3 & 3 & 0 \\
 0 & 3 & 3 & 0 & 0 & 3 & 3 & 0 \\
 0 & 1 & 1 & 0 & 0  & 1 & 1 & 0 \\
 0 & 2 & 2 & 0 & 0 & 2 & 2 & 0 \\
 0 & 0 &  0 & 0 & 0  & 0 & 0 & 0 \\
        \end{matrix}
       \right)
$$

Now instead of one unique number for edge like we had before with $0.5$ we get a gradient in $x$ direction $[1, 1, 2, 2, 3, 3, 2, 2, 1, 1]$.

Now my Questions:

1) Which approach is better, convolving the image with a $\mathbb R^{1 \times 3}$ -Filter or with a $\mathbb R^{3 \times 3}$-Filter?

2) And why is one better than the other?

Thanks

Non-zero pixels after derivation provide you with two informations: a potential location (where the pixel is), and a potential strength (in magnitude) of an edge-pixel. Linear edge detectors are linear filters. A $3\times3$ kernel adds a smoothing effect in the orthogonal direction. It will favor edges with a certain spatial extend, over short-lengthed isolated pixels that do not qualify for real-edges. In other words, an isolated $[0,1,0]^T$, a duplet $[1,1,0]^T$ will be affected a smaller magnitude than three vertical ones $[1,1,1]^T$ after a $3\times3$ kernel.

In different words, the Prewitt  $3\times 3$ separable gradient/smoothing kernel can be written as a (Cartesian) product of two $3\times 1$ and $1\times 3$ kernels:

$$
\left(
\begin{matrix}
 1 & 0 & -1 \\
 1 & 0 & -1  \\
 1 & 0 & -1 \\
\end{matrix}
\right) =
\left(
\begin{matrix}
 1  \\
 1   \\
 1 \\
\end{matrix}
\right) \times \left(
\begin{matrix}
 1 & 0 & -1 \\
\end{matrix}
\right)
$$

This scheme can be generalized: any smoothing filter in one direction, combined with any derivative in the other (orthogonal) direction, can form an horizontal or vertical square/rectangular image edge detector.  You can find a source in Smoothed Differentiation Filters for Images, Meer and Weiss, 1992:

  For the two-dimensional case we restrict ourselves to square
  neighborhoods. In such neighborhoods it is always possible to define a
  separable two-dimensional orthonormal basis built by the Cartesian
  product of two identical one-dimensional bases.

and in Derivative-based Operations or the wiki page Image derivatives
.

If you have clear vertical edges, with same pixel amplitudes, little noise, the $1\times 3$ filter will do the job. Otherwise, good $3\times 3$ are almost always preferable. An image is really a 2D thing, and processing it through 1D tools often wastes information.


How many non-zeros samples in convolution If I have $b[n]$ of length 30 (30 non-zero samples) and $c[n]$ of length 40 (40 non-zero samples).

How many non-zero samples will $a[n]=b[n] * c[n]$ have? (Note '$*$' is a convolution).

I think there will be 30 non-zero samples. Am I correct?

Actually, you can't say how many non-zero values $a[n]$ is going to have because, depending on the values that $b[n]$ and $c[n]$ take, $a[n]$ may equal 0 for some values of $n$ or not.

Nevertheless, I think that you want to know the length of the sequence $a[n]$. As a general rule, for two sequences of length $P$ and $L$, the linear convolution between them returns a sequence of length $N=P+L-1$. In this case, the maximum length of $a[n]$ would be $N=30+40-1=69$ (it can be less than 69 if $b[n]$ and $c[n]$ are such that the convolution between them returns null values for the first or the last values of $a[n]$).

The reason why the resulting sequence has a length of $N=P+L-1$ can be easily seen by doing the convolution graphically. Let's assume for simplicity (although the reasoning without these assumptions would be exactly the same) that the first non-zero value of $b[n]$ appears at $n=0$ and that the same thing happens with $c[n]$. Let's remember what the definition of convolution is:

$$a[n]=b[n]*c[n]=\sum\limits_{m=-\infty}^\infty b[m]c[n-m]$$

This means that we should have in mind $b[m]$ and $c[-m]$, and then start moving $c[-m]$ from $-\infty$ until we get our first non-zero value for $a[n]$. This happens for $n=0$. Then the two sequences always return some information until the last value of $b[m]$ is in the same spot as the first value of $c[n-m]$. This happens for $n=L+P-2=68$. So that leaves us with the length I wrote before: $N=69$.


Complex numbers in a frequency domain of a 2D image I am try to grasp the idea of frequency domain for images. I think get the basics, but now I've stuck with a question that I can't find appropriate answer anywhere.
How are frequency domain and complex numbers in relationship?
I've read this article from DSPguide and I understand that after applying DFT to an image in spatial domain, we get two planes: real and imaginary plane. By doing some calculations, we can get amplitude and phase planes. For what is the phase plane used for?

The real and imaginary plane are also the symmetric and anti-symmetric planes.  Image content that is symmetric around the center lines will show up in the real plane.  Image content that is opposite in brightness across the center lines will appear in the imaginary plane.  Data that is neither or some combination of both will end up in both, with the phase angle indicating the mix or a shift (or rotation wrapping around the edges) from symmetric to anti-symmetric.


Lighting for computer vision automisation Currently I am working on a computer vision project where we want to gather data from a composite material production environment.

I have an example picture here:

What I learned so far that my lighting is not good. It's coming from a fluorescent lamp from that's right above the camera. Apparently I would have to diffuse the light with some sort of material. The reflections of the light are not doing us well.

The red and blue tape was placed there as a test to see if I can segment the image, remove the rest of the image outside of the tape. So that way I can analyse the area inside the tape. I got that working by converting the image to 8 bit and drawing a thick black line and then make a bounding box around that. But with the tape it's not working, presumably due to the reflection. The best outcome would be that the tape is recognized as somewhere between values 0-10.

Basically what I did here:

How can I recreate this in the first image I showed?

What I learned so far that my lighting is not good. It's coming from a fluorescent lamp from that's right above the camera. Apparently I would have to diffuse the light with some sort of material. The reflections of the light are not doing us well.

A "Ring Light" type of light will cast uniform light and minimal shadows to your object. It is used regularly in photography. The link to the product is for demonstration, not to suggest that specific product.

Of course, the other thing with the ring light is that the camera shoots through the ring and it's not casting shadows.

Using diffusers is not efficient for relatively large distances (a few meters), or rather, using a diffuser means that you have to use a much stronger light behind the diffuser to flood a large area.

Now, once the light is flooding the area uniformly, place a sample tape rectangle and extract the Red, Green, Blue (RGB) values of the pixels occupied by the tape. To determine your threshold do a histogram on the RGB channels and try to characterise the modes on the R and the B channel.

A more involved way of doing this is to scan the image once and convert the RGB values to Hue Saturation Value (HSV). You can then do one histogram on the H channel which will have two modes, one on the angle that corresponds to Red and one on the angle that corresponds to Blue. Consequently, you now have one variable (H) with two thresholds and two deviations around the thresholds.

Hope this helps.


Resampling: how many samples to zero-stuff or down sample? I am a little confused about how resampling works with one filter. I get why you choose the lower nyquist cutoff when you combine both interpolation and decimation for resampling by a rational fraction (you don't want to alias and you only get one shot to take care of it when you combine interpolation+decimation).

I will use a practical example: $44.1 \, \textrm{kHz}$ to $48\, \textrm{kHz}$.

$(44100 / 48000) = (441 / 480) = (147 / 160)$

Combining $147$ decimation factor by $160$ interpolation, you can use one anti-alias:
at $44.1/2= 22050\,\textrm{Hz} $ cutoff for a lowpass FIR. However, you still need to zero stuff, right? But how much? 160 zeroes?

Going the other way around $48 \, \textrm{kHz}$ to $44.1\, \textrm{kHz}$, how many samples do you downsample? Certainly not 147 zeroes?

Is $160 - 147 = 13$ zeroes the answer?

In the multi-stage solution, here's how this would work (I think): Stuff $160$ zeroes, anti-alias lowpass at $24\, \textrm{kHz}$. Downsample the output by $147$, antialias lowpass at $22050 \, \textrm{Hz}$. But wait, you're not computing outputs for a lot of those zeroes, so why not just combine the two? Does this leave me with $13$ zeroes for resampling up and $13$ samples to skip over when resampling down before applying the lower of the two cutoffs?

Next thing to wrap my head around: polyphase filterbanks.

EDIT: source of confusion is in last comment to answer, interpolate is upsampling followed by anti-image filter while decimation is anti-alias followed by downsampling. cascading these two obviates the need for two separate filters. if this is the case then would this mean that resampling by a rational factor < 1 entails 2 filters always by the nature of the decimation/interpolation process, while a resample factor > 1 entails that 1 filter may be used instead?

DOUBLE EDIT: thinking about it, you can still use 1 filter, but you must interpolate first, then anti-alias, then decimate.

sorry (definitely not sorry) this is how my brain debugs its malformed thoughts on DSP (:

When you up sample by 160, if you wanted to filter twice, you would actually be applying an "anti-image" filter.  The next filter for decimation would be the "anti-alias" filter.  But in you case, you could just use the anti alias filter once and accomplish the same.  I think perhaps you need to consider what zero stuffing means, to put so many zeros between samples and decimate is to throw away so many samples inbetween.

Btw, you do compute samples for the zeros when you filter.  If you want to understand better, try plotting the spectrum of each stage, and what would happen if only using the 13.

If I understand your question properly, you would effective just upsampling by 13, filter, then downsample.  You'd be back where you started


Difference between Gabor filter implementations I was reading a paper and came across this:

I'd like to implement this in matlab, therefore I initially searched for a premade implementation. I stumbled upon this: http://www.mathworks.com/matlabcentral/fileexchange/5237-2d-gabor-filterver123 .
Though they use a different equation. My question now is: what is the difference between this equation and the one from the paper I read ?

Both equations are related by Euler's formula: http://en.wikipedia.org/wiki/Euler%27s_formula
They are mathematically equivalent, it's well explained on the Wikipedia page.

In a practical point of view, $\cos(x)+i\sin(x)$ allows for easy separation of real and imaginary. On the other hand, $\exp(ix)$ makes it easier to deal with Gaussian since $\exp(x)\cdot\exp(y)=\exp(x+y)$.

Hope this helps!


Convolution and Modulation Theorem with Octave Example As far as I know, convolution theorem stated that
$$
x[n] * y[n] = X(e^{j \omega}) Y(e^{j \omega})
$$
$$
X(e^{j \omega}) * Y(e^{j \omega}) = x[n] y[n]
$$

I tried a simple experiment using GNU Octave which code is provided in the following

x = zeros(1,31);
h = zeros(1,31);
x(14:18) = [20 -30 40 -30 20];
h(14:18) = [1/5 1/5 1/5 1/5 1/5];

So now, we have

x =

 Columns 1 through 20:

    0    0    0    0    0    0    0    0    0    0    0    0    0   20  -30   40  -30   20    0    0

 Columns 21 through 31:

    0    0    0    0    0    0    0    0    0    0    0

h =

 Columns 1 through 10:

   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000

 Columns 11 through 20:

   0.00000   0.00000   0.00000   0.20000   0.20000   0.20000   0.20000   0.20000   0.00000   0.00000

 Columns 21 through 30:

   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000   0.00000

 Column 31:

   0.00000

Next, the following code demonstrates whether $ x[n] * y[n] = X(e^{j \omega}) Y(e^{j \omega}) $

conv_result = round( conv(x,h,shape="same") );
freq_product_result = round( ifft( fft(x) .* fft(h) ) );

Then the results are

conv_result =

 Columns 1 through 25:

   0   0   0   0   0   0   0   0   0   0   0   4  -2   6   0   4   0   6  -2   4   0   0   0   0   0

 Columns 26 through 31:

   0   0   0   0   0   0

freq_product_result =

 Columns 1 through 25:

   0   6  -2   4   0  -0   0  -0   0  -0   0  -0  -0   0  -0  -0   0   0  -0   0  -0   0  -0   0  -0

 Columns 26 through 31:

  -0   4  -2   6   0   4

Indeed both have the $\{4 , -2  , 6 ,  0  , 4 ,  0 ,  6, -2,   4\}$ pattern, but the pattern doesn't appear at the same column index. Is my experiment correct? Should both the result be identical? I think I still don't understand the concept completely.

You are making the standard mistake of computing a linear convolution by performing a circular convolution which is implemented with an inverse DFT of wrong size.

Now assume that your signals $x[n]$ and $h[n]$ are of length $L_x$ and $L_h$ samples each. Then the discrete-convolution of $x[n]$ with $h[n]$ would produce the new sequence $$y[n] =  h[n] \star x[n]$$ whose length is $$L_y = L_x + L_h - 1.$$ This convolution is also referred to as the linear convolution of $x[n]$ and $h[n]$.

According to the well known theorem of DTFT the linear convolution of $x[n]$ and $h[n]$ is related to the multiplication of their DTFTs as: $$ h[n] \star x[n] \longleftrightarrow X(e^{j \omega})H(e^{j \omega})$$ where $H(e^{j \omega})$ and $X(e^{j \omega})$ are the theoretical DTFTs of $x[n]$ and $h[n]$.

The practical verification of this theorem is based on the multiplication of DFTs $X[k]$ and $H[k]$, however, which computes a circular convolution instead of a linear one. So one has to be careful in selecting the correct length for both forward and inverse DFT sizes.

In your implementation you have selected a shorter than necessary length so that there will be aliasing. Apart from specific instances aliasing should be avoid. The following MATLAb/Octave code shows the corrected version.

Lx = 12;           % Length of x[n]
Lh = 15;           % Length of h[n]
Ly = Lx + Lh -1;   % Length of y[n] = x[n]*h[n]

x = randn(1,Lx);   % sequence x[n]
h = randn(1,Lh);   % sequence h[n]

y1 = conv( x, h ); % y[n] computed in time domain by direct linear convolution
y2 = real( ifft( fft(x,Ly).*fft(h,Ly) , Ly) ); % y[n] computed by DFT method