How does an image rotation affect the magnitude of a DFT? Note: This is being applied only to images of b/w text
If we're given some image and plot the DFT spectra as

and then we plot the rotated image's spectra

Is there any way we could recover the angle the image was rotated by based on the magnitudes of the DFT?
The only information I could find on how I would approach this was by using a Hough Transform on the spectral image but that would defeat the purpose of trying to recover the rotation from the DFT.

Is there any way we could recover the angle the image was rotated by based on the magnitudes of the DFT?

Yes.
But some of these alternatives are very similar to the key idea behind the Hough transform.
It very much depends on how the "text image" is produced, but for a typical picture of a page of text (black letters on white background with margins), the "strong" coefficients will be vertical to the direction of the lines of text. The reason for this is not so much the text itself but the regular rectangles that are formed between the lines of letters. These appear as a regular square wave, the Discrete Fourier Transform (DFT) of which forms the "line of strong coefficients".
In general, after getting the spectrum, you would apply a threshold so that you retain only those strong peaks.
Without using the Hough Transform, the most straightforward thing to do here would be to fit a line to the points that make it through the threshold and get a vector that gives you the line's orientation.
The line fitting though would be sensitive to maxima that are not on the same line but still manage to go through the threshold. And this might distort the measurement.
This is very straightforward for the Hough Transform because it "converts" lines to points and due to its summation from different view points is much less sensitive to noise. Detecting line slopes with the Hough Transform is simply a matter of detecting local maxima in the Hough space.
But, this can be slow.
You can reduce the amount of work that the Hough Transform has to do by downsampling your original image in a controlled way. To understand this better, consider this:
Another alternative to the Hough Transform is to produce sums over a wide rectangle that rotates about the centre of the image and select the direction with the highest sum. This "rectangle counting" looks like this:

The gain here comes from summing over less orientations and summing over a band that has some width.
This is equivalent to downsampling the original image by a factor that is equal to the width of the band and reducing the range of angles that the Hough Transform has to produce sums over.
Finally, if we take this to the extreme and assuming that only 1 line is isolated from the thresholding operation, we could get away with just 2 "maxima scan lines".
These now look like this:

After thresholding the DFT spectrum, you are left with a handful of strong "peaks" that (ideally) form a line. Therefore, you could scan the perimeter of the thresholded spectrum for a local maximum. The position of that local maximum corresponds to the slope of the spectrum line.
If the maximum of the horizontal scan line is found in the middle, the line would be at 90 degrees (or, the page would be upright or upside down). If the maximum of the horizontal scan line is found in the beginning, the page would be rotated counter clockwise by 45 degrees. The same holds if the maximum occured at the end of the horizontal scan line which now denotes a 45 degree rotation clockwise.
The same applies for the vertical scan line now because there are some extreme rotations that would bring the line of coefficients entirely off the reach of that horizontal scan line. For those, we have the vertical scanline. And so on.
This is extremely fast because it does not rely on rotations at all, you just scan along the x-axis and y-axis, both of which are simple counters. But it does assume that there is only one strong "grid" in the image that will provide that line of DFT coefficients and it would still be sensitive to noise (or non-typical text images).
Ultimately, this comes down to balancing the accuracy of the assessment to the complexity of producing it.
Hope this helps.


Preloading an IIR Filter Is there a mathematical approach to initializing an IIR filter so that turn-on transients are minimized?
Background: I have a high pass Butterworth design that overloads at the start of a sine wave signal.  If I know the sine turns on at phase =0,i.e. at sine(t) = 0, can I preset the filter registers to minimize the overshoot?

This is a common problem, and has an easy solution. The following assumes you are using a cascade of second-order direct-form-1 sections.
 For a Butterworth highpass filter with a small ratio of Fcutoff to Fs, the numerator coefficients are close to [1 -2 1], and therefore the frequency response of the numerator by itself has a large attenuation at low frequencies. This is compensated by the gain of the recursive section (the poles).
When a low-frequency sinusoidal signal is first applied, the input shift register (corresponding to the numerator) is not filled up until the 3rd clock cycle, and therefore the output of the numerator FIR may be large for the first 2 clock cycles, after which it becomes small. These initial large values get captured in the recursive section, and recirculate for some time until they die out.
So the simple solution is to delay the operation of the filter until the input shift register is filled with the first 3 samples.


Discrete Time Fourier Transform (DTFT) cross correlation property I came across this property of the Discrete Time Fourier Transform (DTFT) and I am having a tough time proving it.

In general, consider two real signals $x[n] \: \& \: y[n]$. If
$$ x[n] \leftrightarrow X(e^{jw}) \\ y[n] \leftrightarrow Y(e^{jw}) $$
Then, $$\boxed{ \sum_{n=-\infty}^{\infty} x[n]y[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(e^{jw})Y(e^{jw})dw} $$

However, when I tried proving this, I am getting an extra conjugate:

$$ x[n]*y[-n] = \sum_{l=-\infty}^{\infty} x[l]y[n+l] \; $$ (Cross Correlation between $x[n]$ and $y[n]$)

$$Also, \; y[-n] \leftrightarrow Y(e^{-jw}).\: As \: x^{*}[n] = x[n], we \: have \; \; Y^{*}(e^{jw}) = Y(e^{-jw})$$

Thus from the convolution property of DTFT, we get:

$$x[n]*y[-n] =  \sum_{l=-\infty}^{\infty} x[l]y[n+l] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(e^{jw})Y(e^{-jw})e^{jwn}dw $$

Setting n = 0 and using the conjugation property as mentioned above (for real signals), we get:
$$\sum_{l=-\infty}^{\infty} x[l]y[l] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(e^{jw})Y^{*}(e^{jw})dw $$

Which is nothing but a result of the cross correlation property of DTFT.

Does the property (mentioned above in the box) hold only for special cases of $x[n]$ and/or $y[n]$ (like the signals being real and even etc.,), or am I making a mistake somewhere?

Thank you.

Your first formula in your question is generally wrong, that's why you can't prove it. The correct formula is

$$\sum_{n=-\infty}^{\infty}x[n]y^*[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\omega})Y^*(e^{j\omega})d\omega\tag{1}$$

which is just Parseval's theorem.

If $x[n]$ and $y[n]$ are real-valued, $(1)$ can be written as

$$\sum_{n=-\infty}^{\infty}x[n]y[n]=\frac{1}{2\pi}\int_{-\pi}^{\pi}X(e^{j\omega})Y(e^{-j\omega})d\omega\tag{2}$$

which is the form you've correctly obtained.

In the special case that $y[n]$ is not only real-valued but also even, i.e., $y[n]=y[-n]$, then $Y(e^{j\omega})=Y(e^{-j\omega})$ holds, and $(2)$ turns into the first formula in your question.


Question about ideal filters Following shows the discrete time Fourier transform of an ideal low pass filter with cutoff frequency $\omega_c$:
$$H\left(e^{j\omega}\right) =
\begin{cases}
1,  & \text{if $|\omega| \le \omega_c$} \\
0, & \text{otherwise}
\end{cases}
$$
Taking inverse discrete time Fourier transform would result in following infinite support impulse response:
$$h[n] = \frac{\omega_c}{\pi} \operatorname{sinc}\left(\frac{\omega_c n}{\pi}\right)$$
If we want to filter out high frequency components of a finite support sequence $x[n]$ in time domain, we have to compute the convolution product of $x[n]$ and $h[n]$. Since $h[n]$ is of infinite length, it is not possible to implement it for example in a computer program and we must use a truncated version of it. But in frequency domain, we can simply multiply $X(e^{j\omega})$ (i.e. discrete time Fourier transform of $x[n]$) with $H(e^{j\omega})$ and then take inverse discrete time Fourier transform to evaluate filtered version of sequence $x[n]$. In brief, my question is why don't researchers and engineers employ frequency response of ideal low pass filters and instead, they have developed a vast body of literature on filter design?
In other words, $h[n]$ is ideal but its frequency response $H(e^{j\omega})$ is practical and realistic.

In other words, $h[n]$ is ideal but its frequency response $H(e^{jω})$ is practical and realistic.

This is not true. To compute $H(e^{j\omega})$, you must evaluate the summation $\sum_{n=-\infty}^{\infty}h[n]e^{-j\omega n}$. A real system can't do this because $h[n]$ is not finite length.
Say you choose to use a rectangular window $w[n]$ (with length $2M+1$) to force $h[n]$ to a finite length and the summation becomes $\sum_{n=-\infty}^{\infty} h[n]w[n]e^{-j\omega n} =\sum_{n=-M}^M h[n] e^{-j\omega n}$. There is still the problem of $\omega$ being a continuous variable. For example, a computer can represent $h[n]w[n]$ in an array, but how is a computer supposed to store $H(e^{j\omega})$ in memory? It should have to store all function values between $-\pi$ and $\pi$ of which there are infinite.


Understanding of Random Process, Random Variable and Probability Density Function I just wanted to confirm my understanding of a Random Process, Random Variable and the its Probability density Function.

Here is the way that I looked a Random Process/Random Variable:

If we consider a sample space $S$ consisting of $n$ outcomes labelled $s_1,s_2,s_3\ldots,s_n$. Suppose if we perform $n$ trials of experiment, we are sure that the probabilities of outcomes $s_1,s_2,\ldots,s_n$ will change over time. So, here the graph $X(s_1,t)$ denotes the changing probability of event s1 over time. graph $X(s_2,t)$ denotes the changing probability of event $s_2$ over time. Similarly we have $n$ functions of time. The whole set of functions $X(s_1,t), X(s_2,t),\ldots,X(s_n,t)$ constitute a Random Process.

Pictorially represented below:

Now when we observe the Random Process at a specific time $t_k$, that is the value at $X(s_1,t_k), X(s_2,t_k),\ldots,X(s_n,t_k)$, if we denote them by $(a_1,a_2,\ldots, a_n)$. Now the mapping between the outcomes $(s_1,s_2,s_3,\ldots,s_n)$ and its probabilities $(a_1,a_2,\ldots,a_n)$ are collectively called as Random Variable. That is we get each outcome and its probability (For example, Probability of outcome $s_1 \rightarrow a_1$ at time $t_k$, probability of outcome $s_1 \rightarrow a_2$ at time $t_k$).

It seems a Random variable $X(t_k)=\{a_1,a_2,a_3,\ldots,a_n\}$ is a collection of outcomes and its probabilities. So, we can get the probability density function of this Random variable by plotting this probability values $(a_1,a_2,\ldots,a_n)$ in a graph. ($X$ axis being possible outcomes $s_1,s_2,s_3,\ldots,s_n$ and $Y$ axis being its probabilities $(a_1,a_2,\ldots,a_n)$. And this probability density function depends upon the time when a random process is observed.

Could you please let me know if the above understanding of a Random Process/Random variable/Probability density function is correct?

when we observe the Random Process at a specific time $t_k$, that is the
  value at $X(s_1,t_k), X(s_2,t_k),\ldots,X(s_n,t_k)$, if we denote them
  by $(a_1,a_2,\ldots, a_n)$. Now the mapping between the outcomes
  $(s_1,s_2,s_3,\ldots,s_n)$ and its probabilities
  $(a_1,a_2,\ldots,a_n)$ are collectively called as Random Variable.

Not true. Random variable is a variable which is ... random where its random characteristic is characterized by its probability density function (PDF). And it is the PDF that is "mapping between the outcomes and its probabilities".

A random process is simply a collection of random variables. $\left\lbrace X(t) \right\rbrace$ is a random process and $X(t_1), ..., X(t_n), ...$ are random variables.

The random variables $\left\lbrace X(t_n) \right\rbrace$ can be identical distributed, i.e. $f_{X(t_1)}(s) = f_{X(t_2)}(s) = \cdots$; or not, i.e. $f_{X(t_1)}(s) \neq f_{X(t_2)}(s) \neq \cdots$ where $f_{X(t_1)}(s)$ denotes the PDF of random variable $X(t_1)$.

Using your notations, $f_{X(t_1)}(s)$ is the probability such that random variable $X(t_1)$ has value $s_1$ of sample space, i.e. $\mathrm{Prob}(X(t_1) = s_1) = f_{X(t_1)}(s_1)$

Update: for the question in comment
Use your notations and if I understand correctly what you mean, for a given time $t=t_k$, we do experiments, and we have outcomes $s_1, \cdots, s_n, \cdots$ which associated probabilities denoted $X(s_1, t_k), \cdots, X(s_n, t_k), \cdots$.

Thus the Probability mass function (or PDF if the support is continous) is $f_{X(t_k)}(s) = \left\lbrace X(s_i, t_k) \textrm{, } s_i \in \left\lbrace s_1, \cdots \right\rbrace \right\rbrace$


Difference between SNR and PSNR I understood that SNR is the ratio of signal power to the noise power. In terms of images, how the original image is affected by the added noise. In PSNR, we take the square of the peak value in the image (in case of an 8 bit image, the peak value is 255) and divide it by the mean square error. The SNR and PSNR are used to measure the quality of an image after the reconstruction. I understand that higher the SNR or PSNR, the reconstruction is good. What I don't understand is how SNR and PSNR differs in terms of their conclusion about the reconstructed image.

What the PSNR of an image concludes that the SNR of the same image can't conclude ?
Simply how the conclusion of PSNR differs from the conclusion of SNR?

Let's start with the mathematical definitions.

Discrete signal power is defined as
$$P_s = \sum_{-\infty}^{\infty}s^2[n] = \left|s[n]\right|^2.$$

We can apply this notion to noise $w$ on top of some signal to calculate $P_w$ in the same way. The signal to noise ratio (SNR) is then simply
$$P_{SNR}=\frac{P_s}{P_w}$$

If we've received a noise corrupted signal $x[n] = s[n]+w[n]$ then we compute the SNR as follows

$$P_{SNR}=\frac{P_s}{P_w} = \frac{P_s}{\left|x[n]-s[n]\right|^2}.$$

Here $\left|x[n]-s[n]\right|^2$ is simply the squared error between original and corrupted signals. Note that if we scaled the definition of power by the number of points in the signal,  this would have been the mean squared error (MSE) but since we're dealing with ratios of powers, the result stays the same.

Let us now interpret this result. This is the ratio of the power of signal to the power of noise. Power is in some sense the squared norm of your signal. It shows how much squared deviation you have from zero on average.

You should also note that we can extend this notion to images by simply summing twice of rows and columns of your image vector, or simply stretching your entire image into a single vector of pixels and apply the one-dimensional definition. You can see that no spacial information is encoded into the definition of power.

Now let's look at peak signal to noise ratio. This definition is

$$P_{PSNR}=\frac{\text{max}(s^2[n])}{\text{MSE}}.$$

If you stare at this for long enough you will realize that this definition is really the same as that of $P_{SNR}$ except that the numerator of the ratio is now the maximum squared intensity of the signal, not the average one. This makes this criterion less strict. You can see that $P_{PSNR} \ge P_{SNR}$ and that they will only be equal to each other if your original clean signal is constant everywhere, and with maximum amplitude. Notice that although the variance of a constant signal is null, its power is not; the level of such constant signal does make a difference in SNR but not in PSNR.

Now, why does this definition make sense? It makes sense because the case of SNR we're looking at how strong the signal is and to how strong the noise is. We assume that there are no special circumstances. In fact, this definition is adapted directly from the physical definition of electrical power. In case of PSNR, we're interested in signal peak because we can be  interested in things like the bandwidth of the signal, or number of bits we need to represent it. This is much more content-specific than pure SNR and can find many reasonable applications, image compression being on of them. Here we're saying that what matters is how well high-intensity regions of the image come through the noise, and we're paying much less attention to how we're performing under low intensity.


Establishing the function of a decaying signal I have a set of samples from a signal. The signal is erratic in behaviour at the beginning but then it settles down and oscillates around the value of 100 (the frequency of oscillation is not necessarily a constant). Having settled down, the signal also begins to decay ever so slightly. Given that I only have samples for the beginning of the signal, how might I estimate the remaining samples in the signal until the point at which it flat-lines. That is, the point at which all consecutive samples will have a value of 100.

Here are my samples:

3
84
200
130
106
80
70
65
66
73
78
92
105
120
120
103
94
79
72
73
78
93
101
109
108
104
90
81
78
79
80
97
101
104
100
92
87
82
87
99

There is Eureqa : Eureqa (pronounced "eureka") is a software tool for detecting equations and hidden mathematical relationships in your data

EDIT : Now below is what I wrote in first place but reading myself again, I haven't really answered the exact question you've asked ... now that I (think) have (above), I just left it because in  my case Eureqa yielded no good results. (pretty much the same issue as you, I was trying to estimate when a turntable would really stop for my DJ-ing software)

So here it is :

What about a http://en.wikipedia.org/wiki/Linear_regression ? You might be interested in the Estimation category, on the right of the page.
There's the http://en.wikipedia.org/wiki/Zero-crossing_rate, it is quite primitive but yielded some nice results for me, you can find out how many times and how far occurences do occur for your value of 100. (You might need a lot more data than a few samples, though ... but you can oversample your signal)
Last thing, you can (and probably should) apply some filtering to the incoming signal, I did use the Moving Average Filter http://www.dspguide.com/ch15/5.htm. Dependings your needs, you can get rid off some of the noise present to even enhance your estimation.

From this great book : The Scientist and Engineer's Guide to
Digital Signal Processing
By Steven W. Smith, Ph.D. It is free ! Chapters of interest for you : 14, 15, 16, 19, 20, 21

Using Wikipedia and this book, I've been able to track a time-code signal for my DJ-ing software.

(if you don't know what a time-code is : http://en.wikipedia.org/wiki/Vinyl_emulation_software)

Also it's a lot of trial and error, magic numbers sometimes, in my case it was 21 samples (for the MA filter) at 44.1Khz, with this I've been able to track changes in the frequency of about 0.25Hz which was way better than when not using the MA filter before the ZCR.


What are some of the disadvantages of parametric PSD estimation methods? Roughly speaking, the methods of power-spectra estimation can be divided in two categories : parametric and non-parametric.

From what I understood, parametric methods assumes that the signal can be more or less modeled (with an AR model, or a behavior + noise model etc..) and try to use this to perform the estimation.

I have noticed that when I use parametric methods, the result I obtain usually looks "nicer" that if I use a non-parametric method, let's say Welch. The thing is, I don't know anything about my signal and do not know if it can be modelized and if it therefore makes sense to trust parametric methods.

This is where my question comes in : What are the know "problems" which are associated with parametric methods ? Can we trust a spectrogram obtained this way even though the signal can't be modelized or doesn't fit the correct framework of the used method ?

Can we trust a spectrogram obtained this way even though the signal can't be modelized or doesn't fit the correct framework of the used method ?

Short answer: anything based on a model will work only when applied to the model. Not quite sure which other answer you were expecting here!

Of course, some "model mismatches" are more benign than others: For example, some classes of estimators are very robust against noise not fitting the original model, but break down when the assumed number of signal components isn't what was specified (example: ESPRIT). Others fail the other way around (Bartlett's method will not like you having negative SNR).


Deriving Frequency Response for 2-pole Zero-Delay Feedback State Variable Filter I have an existing zero-delay feedback (ZDF) 2-pole state variable filter implementation (along the lines of the theory presented in VA Filter Design by V. Zavalishin), and I wish to determine the frequency response so that I can plot it. My usual approach with other filters (e.g. biquads) has been to determine the difference equation by inspection, apply the Z-transform, and then substitute $z=e^{j\omega}$ and plot phase and magnitude against $\omega$ from $0$ to $\pi$.

However I am having difficulty determining the difference equation of this filter.

Here is a diagrammatic representation of the implementation I have:

Here's the implementation of an iteration, corresponding to a single input sample, in pseudo-code. Fs is the sample rate, F is the filter frequency of interest, R is the resonance, x is the input sample, and lp, bp, hp represent the low, band and high-pass outputs respectively:

# z1 and z2 contain state from previous iteration

k = tan(pi / Fs * F)
b = R + k
a = 1 / (1 + b * k)

hp = (x - (z2 + z1 * b)) * a
x1 = hp * k  # temporary
bp = x1 + z1
x2 = bp * k  # temporary
lp = x2 + z2

z1 = x1 + bp
z2 = x2 + lp

For analysis, I assign $v[n]$ to be the output of z1, and $w[n]$ to be the output of z2. From this I am able to determine by inspection the following equations:

$$hp[n] = x[n] - (b.v[n] + w[n])$$
$$v[n] = k.hp[n-1] + bp[n-1]$$
$$w[n] = k.bp[n-1] + lp[n-1]$$
$$bp[n] = k.hp[n] + v[n]$$
$$lp[n] = k.bp[n] + w[n]$$

I'm not sure how to handle such a system of equations algebraically. For example, if I wish to analyse the low-pass response, then I will need to determine a difference equation for $lp[n]$ in terms of $x[n]$. However I have not had success manipulating these equations by rearrangement and substitution even with the help of the Z transform. Is it possible to do this or are there too many unknown variables?

If an algebraic solution is impossible, is there an alternative method to determine the frequency response of this filter?

I do not see any reason why this could not be derived algebraically. But this would be a good application for Mason's gain formula; the references below may be helpful for that approach:

Rick Lyon's Tutorial

Mason's Gain Formula - Wikipedia

Youtube tutorial

Tutorial's Point tutorial


FMCW Radar: Need for Second stage FFT for velocity determination? I have understood that for the purpose of Range determination, the samples of the beat frequency are first assembled per chirp. The beat frequency is directly proportional to the range of the reflecting object.

So if there are N chirps, there will in a 2D matrix be N Rows of M samples each. When the FFT is performed along each row (Chirp), the spectrum analysis yields the value of the beat frequency of the corresponding chirp. Let us assume that the beat frequency is represented by Bin B.

If the object is stationary, then all of the rows will have a significant value and nearly the same value in Bin B [Magnitude and Phase].

If the object is slowly moving away and assuming that each bin represents a large range, then Bin B continues to remain the dominant bin with possibly reducing magnitude and increasing phase.

Question hence is: Why perform a second round of FFT along a column for determining velocity of a range bin?

Why couldn't the individual phase of the FFT samples [Output of first round of FFT] along a column be taken and the velocity deduced from phase difference among bins of the column of a range bin?

After all, the phase of a received signal is proportional to the range of the object.

Actually the second round of FFT is exactly doing that, as by definition differentiation of phase is frequency. In realistic case if target is moving, the target position will not change but its phase will change along the column (also called slow-time, where the row in your case is called fast-time), the second FFT shows this changing phase as Doppler frequency. Furthermore, what second FFT also does it something similar to pulse integration. More the number of rows larger is your integration gain.


Choosing correct filter parameters for IMU sensor datas I have a IMU sensor breakoutboard using ITG3701, LSM303D LinkProduct. This involves accelerometer, gyroscope and magnetometer. I have connected the sensors to a Arduino and sending the datas to another computer over Xbee. My sensor is placed on a wheel along its radius. The controller sends data at every 100 ms delay.

The frequency parameters for sensors , i have chosen to be

// Specify sensor full scale
uint8_t OSR = ADC_8192;  // set pressure amd temperature oversample rate
uint8_t Gscale = GFS_4000DPS; // gyro full scale
uint8_t Godr = GODR_190Hz;  // gyro data sample rate
uint8_t Gbw = GBW_low;  // gyro data bandwidth
uint8_t Ascale = AFS_16G; // accel full scale
uint8_t Aodr = AODR_200Hz;  // accel data sample rate
uint8_t Abw = ABW_50Hz;  // accel data bandwidth
uint8_t Mscale = MFS_2G; // mag full scale
uint8_t Modr = MODR_25Hz; // mag data sample rate
uint8_t Mres = MRES_HighResolution; // magnetometer operation mode

It can be changed to several other allowable limits given here

I want to process the sensor datas through two filters, first HPF (to remove DC components noise) and then LPF (to remove high frequency noise)

I want to do the following

Accelerometer_Raw ->  /HPF/  -> Accelerometer_HPF_datas -> /LPF/ -> Accelerometer_LPF_datas

Similarly

Gyroscope_Raw ->  /HPF/  -> Gyroscope_HPF_datas -> /LPF/ -> Gyroscope_LPF_datas

On the receiver side, i did the following on the sensor datas

   FSa=200; % accelerometer sampling
        FSg=190; %gyroscope sampling

d = fdesign.highpass('N,F3dB',2,3/(FSa/2));
            H_a1 = design(d,'butter');
            H_a1.PersistentMemory=true;
            H_a2 = design(d,'butter');
            H_a2.PersistentMemory=true;
            H_a3 = design(d,'butter');
            H_a3.PersistentMemory=true;
            d = fdesign.lowpass('N,F3dB',2,30/(FSa/2));
            L_a1 = design(d,'butter');
            L_a1.PersistentMemory=true;
            L_a2 = design(d,'butter');
            L_a2.PersistentMemory=true;
            L_a3 = design(d,'butter');
            L_a3.PersistentMemory=true;
            d = fdesign.highpass('N,F3dB',2,3/(FSg/2));
            H_g1 = design(d,'butter');
            H_g1.PersistentMemory=true;
            H_g2 = design(d,'butter');
            H_g2.PersistentMemory=true;
            H_g3 = design(d,'butter');
            H_g3.PersistentMemory=true;
            d = fdesign.lowpass('N,F3dB',2,30/(FSg/2));
            L_g1 = design(d,'butter');
            L_g1.PersistentMemory=true;
            L_g2 = design(d,'butter');
            L_g2.PersistentMemory=true;
            L_g3 = design(d,'butter');
            L_g3.PersistentMemory=true;

Gyroscope_Z axis in degrees per second

Accelerometer_X axis in g

I chosed the cutoff frequncy to be HPF to be 3Hz and for LPF to be 30 Hz. I found there is a very big difference with the gyroscope datas after filtering when the wheel rotates. The output datas are significantly reduced for gyroscope.There are also difference at the peaks of accelerometer datas measured. Hence, when i input this to Kalman, I am getting  lower velocty and the lower angle rotated. I have also verified that the output is wrong, as rotating a complete rotation gives me less than 2Pi radians after using filtering.

My questions

If my method is correct for filtering and choosing the sampling rate, cutoff frequency, order of filter etc.
If not, what filter parameters should I chose instead to get better results?
If HPF for gyroscope is necessary?

Additional:

After suggested the lower frequency limit should be much lower: I did the following High pass filtering at 0.05 Hz, The accelerometer datas are good, however i am somehow filtering out the content from Gyroscope datas

The lowpass filter at 3 HZ gave me much better result.

PS: I found the these 3 and 0.05 HZ from the FFT analysis

1) HPF I would say you probably can't have one.
The  accelerometer has similiar considerations.

The gyroscope stream is in terms of xxx output=yyy angular velocity.  Any HPF analog or digital would sag to zero when you attach it to the hub of bicycle wheel or hub of a motor.

The alternative to a HPF is "sample and hold" a value when the unit is at rest.  There are three alternatives.

1-a) During manufacturing take a string of readings at rest and store the value as an offset.

1-b) During power on do the same thing and tell the operator not to disturb it until a ready signal/light comes on.

1-c) Just live with the error of about 1/2 deg/sec ; or 1 rpm after 12 minutes.   That's as I read the error from your data.

If you are willing to live with the sag then I can do the calculations; but I keep thinking of the motor, or bicyclist, going faster and faster trying to keep up to the "speed".

An "alternative" is to differentiate the input and then integrate but that usually doesn't work out very well due to errors introduced during differentiation.

2) LPF

I got good results with a cutoff frequency .34 hertz and a fourth order bessel filter.  My octave 4 (a matlab like program) program and data are at:

https://www.dropbox.com/sh/av5w0qi4myc2tbm/AACxJllMNDsxT2x9nOvjFjMPa?dl=0
Here are the graphs for the first column of gyros.

)2-a Overlay of filtered output versus raw.

)2-b  Filtered signal

)2-c  Normal probability plot pplot()  Straight lines indicate that the errors are gaussian in distribution.  The lapse in the center is a deadband due to resolution

There are other tests but these exemplify the process.  I always check fittings with pplot (or relatives) to see if the fitting leaves gaussian in distribution.  The "staightness" fits into our visual system and is easy.  If not I want to know why; if so I want to know why :) The point is that for some (most?) process trying to dig into gaussian residues is a thankless task.  The truth is I would also look at the autocorrelation as well; but I doubt if you are interested.  One other thing is that one can also plot expected bounds on "straightness" if you want.  In fact this graph indicates that perhaps the filtering is excessive.


Why there is a peak at zero frequency? I have a very simple code for ploting frequency spectrum of a signal in MATLAB. My center frequency is at 5e6 but I see a much higher peak at zero frequency despite the fact that I remove DC component from my signal. Why this peak happens?

Fs = 40e6 ;
[freq,z_freq] = freq_spectrum(RF_ch,Fs) ;
figure;
plot(freq,abs(z_freq));xlabel('Frequency');title('Frequency Spectrum')

function [f,z] = freq_spectrum(x,Fs)
N = size(x,1) ; %size(x)
f = -Fs/2:Fs/(N-1):Fs/2;
x = detrend(x) ; % or %x = x - mean(x);
z = fftshift(fft(x))./numel(x) ;
end

The OP may have removed DC (the mean of the signal) and that would be clear by zooming in on DC to visibly see the null that does exist there. If the null does not exist, then that would mean the mean signal was not properly removed.

The width (effective bandwidth) of the null if done by a simple subtraction of the mean would be approximately $1/T$ in Hz where $T$ is the duration of the data capture in seconds. This appears to be a very long capture and therefore the null would be very tight.

The plot suggests that there is a lot more low frequency content besides a DC offset. If the expected signal is at 5e6, then a high pass filter with a cutoff at approximately 1e6 would be easy to implement and more effective in this case. Further a log scale plot of the vertical axis will make for a better presentation of the spectrum details.


Applying Window Function to Low Frequency Signals I'm computing 1 sec. hanning windows from a signal that has a variable frequency between 1-4Hz. My first impression is that data loss caused by the distortion due to the window function is greater than the correction that it is expected to give while taking FFT of that window. Am I right?

Is it a bad approach to apply window functions to such windows with size of at most 4 wavelengths?

A 1s Hann window will have (in the frequency domain) a main lobe width of 4Hz. This will severely mask your signal components.

Even a rectangular window will be ~2Hz main lobe.

You should use a greater windowing time depending on the accuracy you need for frequency detection.


Sampling of a band-limited analog signal: DFT and aliasing I saw this question on one of the sites related to DFT:

  The analog signal $x(t)$ is band-limited to $40\textrm{ Hz}$. Suppose the signal is sampled at the rate of $100\textrm{ samples/second}$ and that at this rate $200\textrm{ samples}$ are collected. Then $200$ zeros are appended to the $200\textrm{ samples}$ to form a $400$-point vector. Then the $400$-point DFT of this vector is computed to get
  $X[k]$ for $0 \leq k \leq 399$.

  Which DFT coefficients are free of aliasing?
  The DFT coefficient $X[50]$ represents the spectrum of the analog signal at what frequency $f$?
  (Give your answer in $\textrm{Hz}$).

The answer's provided were:

  All of the DFT coefficients are free of aliasing. The sampling rate is more that twice the maximum signal frequency.
  The DFT bin width is $100/400$ or $0.25\textrm{ Hz}$. The $50^{\rm th}$ DFT coefficient corresponds to the frequency $50$ times
  $0.25\textrm{ Hz}$ or $12.5\textrm{ Hz}$

I'm not really sure if this is right because when you append zeros you are upsampling and the sampling rate becomes $400\textrm{ Hz}$. So when you take a $400$-point DFT the bins are seperated by $1\textrm{ Hz}$.

Please let me know if I am right or not. If I am wrong please explain the answer to me.

Jose.

The provided answers are correct. Appending zeros is not upsampling and therefore the sampling rate is still $R_\mathrm{s}=100\, \text{Samples/s}$. Accordingly, the frequency spacing after taking the N-DFT is $f_\mathrm{d}=(R_\mathrm{s}/N)=0.25\, \text{Hz}$.

Upsampling can be achieved by inserting $N_\mathrm{U}$ zeros after every sample followd by lowpass filtering.


Interpolation of missing audio signal in a video sequence Suppose there is a video sequence and there are some frames for which the audio data is missing. I want to interpolate the missing audio data on the basis of the correlation of the audio signal with the pixels of the frames. Is there any existing research on how to do this effectively?

I was thinking of using the Wiener-Hopf equations, but I can't determine if there is any scalar parameter of the video which is correlated with the audio.

Some insight on this would help.

In general, no.  It's possible to create an MP4 file where the video frames and the audio are completely uncorrelated (edited, recorded from different source events, synthesized using different random seeds, etc.)

However, in some cases, it might be possible to infer a correlation, such as video of a human speaker whose face is clearly visible enough to allow lip reading algorithms estimate some percentage of the phonemes, with a statistically significant match against a portion of the nun-corrupted audio, etc.


Digital Phase Modulation baseband or passband? Say I want to transmit messages using PSK:
I send the message
$$s_m(t) = \text{Re} \left [g(t) e^{\frac{j 2 \pi (m-1)}{M}}e^{j 2 \pi f_c t} \right]$$
where $f_c$ is the carrier frequency. Now my question is,

is this a baseband representation or a passband?
So is $f_c$ the center frequency of a certain channel? Or would we need to modulate $s_m(t)$ to the channel frequency?
So then if we sent $s_m(t) \text{cos}(2 \pi f_b)$ where $f_b$ is the center frequency of the channel we are communicating on, at the receiver we would need to convert the modulated signal back to its low pass equivalent $s_m(t)$ so that we can apply a detection algorithm to demodulate the signal?

Assuming $g(t)$ is a baseband signal, $s_m(t)$ is passband. You can see that it is by multiplying the two complex exponentials, expanding the result into sine and cosine terms, and then taking the real part.
An interesting thing about PSK signals is that they are an example of quadrature modulation, and in consequence, you can represent them either as real passband signals, or complex baseband signals. If you assume $s_m(t)$ is a PSK signal, and given that it is real, it must also be passband.


Find B given A and impulse response of B/A Here is a small MATLAB function that computes B in B/A given A and the impulse response h of B/A

function b = numf(h,a,nb)
%NUMF   Find numerator B given impulse-response h of B/A and denominator A
%   NB is the numerator order.

nh = max(size(h));
impr = filter(1,a,[1 zeros(1,nh-1)]);

egqs = toeplitz(impr,[1 zeros(1,nb)]);

b = h/egqs';

Could someone help me understand why these steps give B.

I get what filter() and toeplitz() do, but not the last line : b = h/egqs'

Thanks for any help,

Jeff

The vector impr holds the impulse response of length nh of an all-pole filter with denominator coefficients given by a. Let's call this impulse response $g(n)$. The impulse response $h(n)$ must satisfy the equation

$$h(n)=(g*b)(n)\tag{1}$$

where * denotes convolution. With the given data, this results in an overdetermined system of linear equations, the system matrix of which is the Toeplitz matrix egqs. The last line in the code solves this overdetermined system of linear equations corresponding to (1) in a least-squares sense. You could achieve the same using the backslash operator:

b = egqs\h';


Methods of scaling FFT power I wanted to know what are the well-known methods of scaling the FFT power of a signal? I am referring to EEG signal but this should apply to all signals I guess.

I need to plot the powerspectrum/PSD of frequencies but without correct scaling lower frequencies are prominent while higher frequencies are not.

What are the math formulas? (i.e.: 10log10(power))?

In displaying the PSD of signals it's common practice the subtract the DC component first;
$$ x_{ac} = x_{tot} - x_{dc} $$

So that spectral leakage will not suppress low magnitude high frequency components.

Furthermore, in some spectral analysis applications one may need a proper scaling between what is computed by the DFT of the windowed sampled sequence $v[n] = w[n]x_d[n]$ and the true continuous time CTFT spectrum $X(\Omega)$ of the analog signal $x_a(t)$.

That scale will be a function of the window length $M$ and sampling frequency $F_s$. Lets see the implications of this from an example:

Consider a continuous-time signal $x_a(t) =  \left( \frac{\sin(W t)}{\pi t} \right)^2$ whose CTFT $X_a(\Omega)$ is
$$ X_a(\Omega) =
\begin{cases}
A (1 + \frac{\Omega}{2W}) &, \text{for} -2W < \Omega < 0 \\
A (1 - \frac{\Omega}{2W}) &, \text{for} -0 < \Omega < 2W \\
\end{cases}
$$

Which is a bandlimited triangular pulse with $A = X_a(0) = \frac{W}{\pi} $ These two signals are plotted below from their formulas for $W = 7\pi$ rad/second.

Now, this analog signal is (critically) sampled with a sampling period of $T_s = \frac{2\pi}{4 W} =\frac{\pi}{2W} $ such that we obtain the sampled discrete-time signal $x_d[n] = x_a(nT_s)$ as
$$ x_d[n] =  \left( \frac{\sin(W n T_s)}{\pi nTs} \right)^2 = \left( F_s \frac{\sin(w_c n)}{\pi n} \right)^2 $$

Where $w_c = W T_s $ and $F_s = 1/T_s$. In particular for the selected critical sampling rate $w_c = W \frac{\pi}{2W} = \frac{\pi}{2}$ and $F_s = \frac{2W}{\pi}$.  At this stage we can consider the DTFT $X_d(e^{jw})$ of the sampled sequence $x_d[n]$ either from its formula as given above or from the relation between $X_d(e^{jw})$ and $X_a(\Omega)$ from the sampling operation as given by:
$$ X_d(e^{jw}) = \frac{1}{Ts} \sum_k X_a( \frac{w-2\pi k}{T_s}) $$

It can be seen (by either methods) that $X_d(e^{j0}) = \frac{1}{T_s} X_a(0) = \frac{W}{\pi T_s}$. With $T_s = \frac{\pi}{2W}$ it becomes $X_d(e^{j0}) = \frac{1}{T_s} X_a(0) = 2(\frac{W}{\pi})^2$ which is $98$ for $W=7 \pi$.

Plotting $x_d[n]$ and $X_d(e^{jw})$ yields the following:

It can be seen that the discrete-time frequecy spectrum is extending from $w=-\pi$ to $w=\pi$ as it was critically sampled.

Now we shall return to obtaning a finite-length set of samples $v[n]$ from $x_d[n]$ by a sutiable window of length $M$ :
$$ v[n] = w[n] x_d[n] $$

From which the DTFT $V(e^{jw})$ is deduced as
$$V(e^{jw}) = \frac{1}{2\pi} V(e^{jw}) \star X_d(e^{jw}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} W(e^{j\theta}) W(e^{w-j\theta}) d\theta$$

At this point we have to specify a particular window type such as rectangular, Hamming or Kaiser etc. Lets choose the simplest such as rectengular window of length $M$ whose DTFT is:
$$W(e^{jw}) = e^{j w \frac{M-1}{2}} \frac{ \sin( w M/2) }{ \sin(w/2) }$$
Note that the peak value of $|W(e^{jw})|$ is $M$

Lets evaluate the result of this convolution at the frequency $w=0$:

$$V(e^{j0}) = \frac{1}{2\pi} \int_{-\pi}^{\pi} W(e^{j\theta}) W(e^{-j\theta}) d\theta$$

The exact integration is difficult. But proivded that the window length $M$ is long enough so that the spectrum of the window $W(e^{jw})$ is highly concentrated in its main lobe which is between $-\frac{2\pi}{M} < w < \frac{2\pi}{M}$ and has a width of $\Delta_{ML} = \frac{4\pi}{M}$ , we can get an approximation for $V(e^{j0})$ as

$$ V(e^{j0}) \approx 0.5 \frac{1}{2\pi} X_d(e^{j0}) W(e^{j0}) \frac{4\pi}{M} = X_d(e^{j0})$$The approximation gets better as $M$ gets larger.

Hence we conclude that the scaling for the DTFT $V(e^{jw})$ is the same as that of $X_d(e^{jw})$ for a rectangular window of long enough length $M$ (the scale may be modified with the sum of window length for other windows lengths)

Finally we shall consider the DFT $V[k]$ associated with the DTFT $V(e^{jw})$.
Since for finite length signals (such as $v[n] = w[n]x_d[n]$ of length $M$) the DFT values $V[k]$ are uniform samples of the associated DTFT $V(e^{jw})$ there is no further scaling due to DFT stage. We conclude that:

$$\boxed{ V[0] = V(e^{j0}) = X_d(e^{j0}) = \frac{1}{T_s} X_a(0) = \frac{W}{\pi T_s} }$$


Does any paper use energy harvest to in the baseband? Most of RF energy harvest system just harvest $\rho$  of the RF signal for charging the battery,and the other part RF energy will be processed to baseband signal by RF Chain to do the information decoding,like this below

Now i have a new thinking,if we receive signal and use the RF chain to let the signal become the baseband signal first,then use the power splitter to let $\rho$ of baseband signal power to charge the battery(energy harvest).Can we do that?

Is there any paper which use this method to do the energy harvest?like this below

No.

Your question neglects the fact that a mixer can't transfer the whole received power to baseband (mixers are typically active, and most of the baseband energy comes from the LO that the device has to generate itself).

You can use any nonlinear device, even one without an external power source, but the mixed output's energy will always be way lower than the energy of the inputs, simply by the pure math behind what mixing does.

In general, when asking such question, you'd start with "where does what energy come from?", and with a clean, and simplified, block diagram of the receive RF chain including the mixer, and that would have answered that. If you're researching power harvesting, this is the MAIN thing you'll need to understand!


Amplitude of the fft is less than what is expected This is my MATLAB code
Amp = 0.3;
freqHz = 1000;
fsHz = 160000;
dt = 1/fsHz;
t = 0:dt:2000*dt;
sine = Amp * sin(2*pi*freqHz*t);
t = 0:2000
subplot(2,2,1)
plot(t,sine)

y1 = sine;
Fs = 160000;                   % samples per second
N = length(y1);            % samples
dF = Fs/N;                 % hertz per sample
f = -Fs/2:dF:Fs/2-dF + (dF/2)*mod(N,2);      % hertz
Y1_fft = fftshift(fft(y1))/N;
subplot(2,2,2)
plot(f , (abs(Y1_fft)) - min((abs(Y1_fft))));
legend(num2str(i))
shg;

and this is the result:

The amplitude of my sin is 0.3 so I expect to get the FFT signal with the amplitude of 0.15 at 1Khz while I get 0.096 as the amplitude. My goal is to extract the amplitude of a tone among different frequencies.

The amplitude peak is hidden between FFT result bins for frequencies that are not exactly integer periodic in the FFT's length.  But you can interpolate to find the results for those frequencies.  Try windowed Sinc interpolation of the real and imaginary components, then compute the magnitude of the complex result.  If you don't know the exact frequency, you can use successive approximation to find an interpolated magnitude peak between FFT result bins.


How do I apply the Nyquist criterion on the following spectrum to see if it satisfies the criterion or not? I am trying to find out if the pulse in the figure satisfies the Nyquist criterion for a modulation interval T. I'm used to the Nyqusit frequency, but I cannot grasp what the Nyquist criterion is about.

The Nyquist criterion for avoiding intersymbol interference (ISI) requires the pulse to have zeros at multiples of the symbol interval $T$, such that symbols do not interfere with each other. The time domain requirement is

$$h(t)=\begin{cases}1,&t=0\\0,&t=nT\end{cases}\tag{1}$$

In the frequency domain, criterion $(1)$ can be formulated as

$$\sum_{k=-\infty}^{\infty}H\left(f-\frac{k}{T}\right)=T\tag{2}$$

where $H(f)$ is the Fourier transform of $h(t)$. I.e., the shifted spectra (shifted by the symbol rate) must add up to a constant.

I leave it up to you to decide whether or not the given pulse satisfies $(1)$ and $(2)$.

I suggest you go through your lectures notes because it would be surprising if the Nyquist criterion has never been mentioned and you get an exercise like this.


How many things can they be done with a spectrogram? I found some things that can be done with a spectrogram.

filter frequencies by setting the bins to zero
observe what frequencies make up the signal.
observe the energy or
amplitude of each frequency, the whitest pixel has more volume and the black pixels have little amplitude.
observe harmonics

It does not occur to me that other things can be done.

A spectrogram is a reversible (and generally redundant and complex) transformation that is used to separate components, unmix latent signals, some of which being related to frequency analysis (Nota: sometimes, this name is used for the (squared) magnitude of the above. I do use it as a short-hand for the breeds of short-term Fourier transforms).
We can note that the classical Fourier transform is a (simplistic)   form of spectrograms. So anything you can do with Fourier, you can expect to do it (possibly better) with spectrograms.
For instance, spectrograms can be used  to sparsify representations (reduce the proportion of useful components) or to emphasize features from non-stationary signals. So any signal processing technique, that could use the former sparsifying preprocessing (could be an FFT alone)), fits the picture: time-frequency cancellation, detection, noise estimation and removal, segmentation, deconvolution, adaptive filtering, source separation, learning, etc.
For instance, even mp3 or JPEG compression are based on a sort of (non-redundant) spectrograms.


DSP processor: some beginner's questions I have started making a project involving FFT to some accoustic signal(so far I completed some theory), with a DSP processor, and according to some calculations, make the output going to a bluetooth (or any other wireless device) chip, which will transmit to a bluetooth receiver.

So far I have seen wireless chips, and DSP proccessors, but not from the same company.

Programming each chip that came from a different company requires different softwares or big differences in commands?
Is there a company making both DSPs and bluetooth?
If you have some tips regarding the chips I will be happy to read them :) (easiest way to begin, where I should go into the architecture and what is less importent and can look at it from a bigger picture)

I'v have some experience with TI C55 and C6000 series DSPs, I'll give you some points, hope to be helpful.

  Programming each chip that came from a different company requires different softwares or big differences in commands?

Absolutely! Although the code body might be the same, but since every company has its special processor architecture, their compilers would be different. Even two DSPs of a same company might be too different, for example Texas Instruments C5000 processors are too different from C6000, since their memory architecture is different. I mean, you can't simply use the same segment of code for both.

  Is there a company making both DSPs and Bluetooth?

Usually most of the DSPs are only Microprocessor and a few essential Interfacing (like I2C, SPI,...) are considered for them, however some are more sophisticated System on a Chip designs might include Bluetooth. Basically you can buy and add a Bluetooth module to any DSP you want.
https://www.sparkfun.com/products/12574

  If you have some tips regarding the chips I will be happy to read them :) (easiest way to begin, where I should go into the architecture and what is less importent and can look at it from a bigger picture)

It depends on your programming level and those calculations you've mentioned above, however, I think the best chip to begin with is MicroChip DSPIC, which is simply a PIC micro associated with some DSP blocks (like circular conv, FFT ,etc), I think your last choice should be TI's DSPs or AD's Shark DSPs, they are too complicated and frustrating to begin with.


Fourier Transform frequency translation property in practical applications From the material of my Signal Analysis and Processing course:

  The frequency translation property is extremely important for
  practical applications. It allows taking a signal from a portion
  of the spectrum and relocating it onto any other portion of the
  spectrum. Thanks to this property, for instance, radio
  transmission is possible, including FM radio, cell phones, TV
  transmission and so on.

I understand this property but I can't find anything on the Internet regarding its practical applications. Do you have links to any resource of this type?

I can't find anything on the Internet regarding its practical applications.

Translation, in general, occurs when a series is multiplied with (or "modulated by") a sinusoid. I write "series" because the translation effect can be applied to both the time domain and the frequency domain.

The main thing to keep in mind here is that when you multiply a time series with a sinusoid you "shape" its spectrum in a predictable way.

This is the standard trick used by Heterodyning which mentions a few more applications of the same concept at its "Applications" section.

At the heart of this property are the Product To Sum (and back) trigonometric identities. These point back to the subject of modulation and phase detection too.

But also, this property, along with a similar property that holds for logarithms were behind the proliferation of various types of slide rules which were early aids to support quick calculation. Being able to reduce a complicated calculation involving trigonometric functions down to a set of much simpler additions fuelled rapid expansion in a number of different disciplines in the past. Slide rules used to be found in all sorts of ships.

Hope this helps.

EDIT:

Another variation on the same theme is the modification of the low pass digital filter.

The simplest digital filter everyone is introduced to in DSP is the Sinc Filter. A practical sinc filter, is a low pass filter. By modulating a sinc filter with different types of sinusoids (real or imaginary) your can change its cut-off frequency (in real time) or transform it to any other type of filter (That is, from low pass you can transform it to a high-pass, band-pass or band-reject. One set of coefficients, different modulations).


How to create matched "minimum phase" for a system of parallel FIRs? Problem statement

I have a collection of magnitude (only) responses I'd like to turn into FIR filter kernels that are

matched in phase
"minimal" in phase, with respect to the complete system
have little to no "pre-ring"
maintain the design magnitude to a high degree

In the ideal case, I'm hoping for a solution that is analytic rather than parametric.

Meeting the first requirement can be expressed as:

$$
\phi_s(f) = arg\{ H_0(f) \} = arg\{ H_1(f) \} = ... = arg\{ H_n(f) \}
$$

And, is very easy to meet in a linear phase case.

Accumulating phase

In attempting to meet the second condition, I've tried something akin to a "cascading" of phase given minimum phase versions of each filter. I've done this by

returning the "Hilbert" minimum phase for each magnitude response
accumulating these responses, and then defining this as the system phase
reassigning system phase to each response

Here's an expression for the accumulation I used to find the system phase:

$$
\phi_s(f) = \sum_{i=0}^n - \mathscr H \{ \ln( |H_n(f)| ) \}
$$

Unfortunately, for my responses, this results in a poor outcome, returning kernels with way too much smearing / lag at low frequencies in the time domain. All this smearing is far from "minimal" in any sense!

I'm suspicious this may be a result of the two definitions of minimum phase discussed earlier.

Minimum phase of L-2 norm of magnitudes

The most promising results I've had after a bit of hacking is to take the L-2 norm of my magnitude set and use this as a minimum phase target:

$$
|H_{L2}(f)| = \biggl(\sum_{i=0}^n |H_n(f)|^2 \biggr)^{1/2}
$$

$$
\phi_{L2}(f) = - \mathscr H \{ \ln (|H_{L2}(f)|) \}
$$

Setting $\phi_s(f) = \phi_{L2}(f)$ returns kernels with coefficients packed towards the beginning and end of the table. Rotating the table, a linear phase offset, by $\max\bigl(\phi_{L2}(f)\bigr)$, turns the end into a "pre-ring".

While not beautiful, this does get me closer to what I think I'm looking for.

Some more details on my magnitudes

$|H_0(f)|$ is allpass; the remaining are highpass with every increasing roll-off slopes and cut-off frequencies.

Questions

Is there a nice off the shelf know solution? In my mind, the problem is very similar to a matched phase cross-over.

Am I asking too much? In the ideal case I'd like to have no "pre-ring". Maybe that's not possible?

Am I asking the right question(s)? In other words, is "minimum phase" what I'm looking for. What I really want is for my kernels to be matched in phase, and respond "as soon as possible".

The Hilbert "transform" or relation define one phase response for a given magnitude response, so you can't get both matched and minimum phase in your case.

Regards.


Struggling to find a good definiton for Gabor-Descriptor I am attempting to implement the Gabor Descriptor based on the following paper Image Copy Detection Using a Robust Gabor Texture Descriptor. It is the best paper that I have found on the subject. It is not perfect, as it has some ordering issues with some equations and paragraphs that only make sense after later sections. No biggie, it is only a matter of time consumed to decipher the intent.

What is stopping me right now, is that one of the parameters is never defined in the paper. It first appears in $(1)$ of Section 2. There we have the following factor: $e^{2\pi jWx}$

$W$ is later assumed to be 1 and $x$ is the parameter of the function, but $j$ is never defined (though it is used again in several subsequent equations, especially in section 3.1, of most interest is $(11)$).

So I am seeking either a definition of $j$. A nice bonus would be a better paper (or book) defining the Gabor Descriptor.

In this context, $j$ is the imaginary unit: $j = \sqrt{-1}$. The complex exponential function represents what they call an "oriented sinusoidal grating." The fact that $j$ is the imaginary unit in this case can be verified by noticing the relationship between $(1)$ and $(2)$ in Section 2; this is an example of the frequency-shifting property of the Fourier transform. That is, multiplying the Gaussian function $h(x,y)$ by the complex exponential function results in a shift in the Fourier transform of the result along the axis that the sinusoid runs across (the $x$ axis in your example).


Undesirable zeros in sampled signal using train of impulses I'm trying to analyze a signal like the example below, but it contains undesirable zero value samples, is there any way to eliminate them? I'm using Matlab.

Actually, the problem is this: in a sampling process I am using a train of impulses but that train I formed it with the cardinal sine function and it has generated values with zeros. Then, that train multiplies it in time to the sinusoidal signal and obtain samples (samples with zeros).

k=5;
t=0:Ts:T;
for n=min(t):k*Ts:max(t)
    pdirac=pdirac+sinc(100*(t-n));
end
signal=signal1.*(pdirac);

If you just dont want to see them in your figure, just replace them with nans. Lets say that your signal is s. Then

s(s==0)=nan;
plot(s);


Synthesis discrete time signal from fourier coefficients Following information is given about a signal $x[n]$

$x[n]$ is real and even signal
$x[n]$ has a period $N=10$ and Fourier coefficients $a_k$
$a_{11}=5$
$\frac1 {10}\sum_{n=0}^9 |x[n]|^2=50$

How can we obtain $x[n]$ from these information?
I know the first info makes the Fourier coefficients real and even, third info makes the Fourier coefficient $a_1=a_{11}=5$(since Fourier coeffs are also periodic with N=10) and fourth info is Parsevals theorem.

As you indicated, Fourier coefficients are even from first point. This tells you that $a_9 = a_1$, so $a_9 = 5$.

Then from Parseval's theorem, you get

$$
   \begin{align}
    \frac{1}{10} \sum_{n=0}^9 \left|x[n]\right|^2 &= \sum_{k=0}^9 \left|a_k\right|^2 \\
     &= 50
   \end{align}
   $$

Since $a_1^2 + a_9^2 = 50$, the other coefficients must be 0.

Then $x[n]$ can be computed as:

$$
   \begin{align}
    x[n] &= \sum_{k=0}^9 a_k \exp(2\pi j n k / 10) \\
         &= a_1 \exp(2\pi j n / 10) + a_9 \exp(2\pi j n \cdot 9 /10) \\
         &= a_1 \exp(2\pi j n / 10) + a_9 \exp(-2\pi j n /10) \\
         &= a_1\left[\exp(2\pi j n / 10) + \exp(-2\pi j n /10)\right] \\
         &= 2 a_1 \cos(2\pi n/10) \\
         &= 10 \cos(2\pi n/10)
   \end{align}
   $$


Independence of Noise at Each DFT Output My math may be a little rusty, so I would like confirmation or correction or my calculations here.

Given white noise samples, $x_i$, which are IID and zero-mean, and variance $\sigma^2_x$. I want to show that the outputs of any pair of DFT filters is also IID.

Knowing covariance of the outputs should be enough to draw a conclusion. So I examine the following:

Given the outputs $X_k$ and $X_{k'}$ defined as

$$X_k=\Sigma_{n=0}^{N-1}x_ne^{-2\pi ikn/N}$$

$$X_{k'}=\Sigma_{m=0}^{N-1}x_me^{-2\pi i{k'}m/N}$$

We can calculate the expected value of $X_k X^*_{k'}$

$$\begin{eqnarray*}
E(X_k X^*_{k'})
&=&E\left(\sum_{n=0}^{N-1}x_ne^{-2\pi ikn/N}\sum_{m=0}^{N-1}x_me^{2\pi i{k'}m/N}\right)\\
&=&E\left(\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}x_n x^*_m e^{-2\pi i(kn-{k'}m)/N}\right)\\
&=&\sum_{n=0}^{N-1}\sum_{m=0}^{N-1}E\left(x_n x^*_m \right) e^{-2\pi i(kn-{k'}m)/N}\\
\end{eqnarray*}$$

We know that $E(x_k x^*_{k'})=\sigma^2_x$ when $k={k'}$ and $E(x_k x^*_{k'})=0$ when $k\neq {k'}$. So we can simplify:

$$\begin{eqnarray*}
E(X_k X^*_{k'})
&=&\sigma^2_x\sum_{n=0}^{N-1}e^{-2\pi i(k-{k'})n/N}\\
&=&N\sigma^2_x|_{k={k'}}
\end{eqnarray*}$$

But if I'm right so far, I still am fuzzy on the case where $k \neq {k'}$. I expect it to equate to zero. Can someone help me out with that part?

Let $m=k-k'$. Then the sum in your final equation becomes

$$\sum_{n=0}^{N-1}e^{-2\pi imn/N}=\frac{1-e^{-2\pi im}}{1-e^{-2\pi im/N}},\quad m\neq 0\tag{1}$$

where I've used the formula for the geometric series. Noting that $e^{-2\pi im}=1$ for any integer $m$ shows that the sum in $(1)$ equals zero.

[Note that your sums should only have $N$ elements, i.e., their indices should go from $0$ to $N-1$ (instead of $N$).]


A (Clear) Line Between Computer Vision and Image Processing I have been working in and learning about both Computer vision and Image processing for a few years now, and I believe that I am not a complete beginner any more.

Still, after all these years, it is hard for me to tell for any particular part of my work whether it is mostly Computer vision related, or if it is Image processing. I just can't see the line -- when I work, study and research, I read reference materials with both keywords.

So, I am interested in the definition of the Computer vision and Image processing fields, with focus on the comparison (differences as well as overlaps) of the fields.

Additionally, I think it would be valuable to have an example of a (conceptual or existing) practical applications, projects and works dealing with/utilizing:

solely (or mostly) Computer vision tools and ideas
solely (or mostly) Image processing tools and ideas
a combination of tools and ideas from both fields

with special attention paid to what makes it one and not the other, or what makes it both.

I understand that these fields are highly related and that the "line" might not be as clear as this question is asking for, but I hope you understand that the point of this question was not to design a simple decision rule for classifying (my) work, but rather a better understanding of the focus and goals of these fields. Also, any additional information that seems on-topic with my question tone is welcome, even if I did not specifically ask for it.

I believe Gonzalez and Woods are competent enough to rely on their opinion:

  There is no general agreement among authors regarding where image processing stops and other related areas, such as image analysis and computer vision, start. Sometimes a distinction is made by defining image processing as a
  discipline in which both the input and output of a process are images.We believe this to be a limiting and somewhat artificial boundary. For example, under this definition, even the trivial task of computing the average intensity of an image (which yields a single number) would not be considered an image processing operation. On the other hand, there are fields such as computer vision whose ultimate goal is to use computers to emulate human vision, including learning and being able to make inferences and take actions based on visual inputs. This area itself is a branch of artificial intelligence (AI) whose objective is to emulate human intelligence. The field of AI is in its earliest stages of infancy in terms of development, with progress having been much slower than originally anticipated. The area of image analysis (also called image understanding) is in between image processing and computer vision.

So I would say, that primary difference is in goals, not methods. For example, if the goal is to enhance image for later use by humans, than this may be called image processing. And if the goal is to emulate human vision (be it object recognition, defect detection or automatic driving), then it is closer to computer vision. Note, however, that emulating human vision by definition may also require image enhancement, so in most real cases computer vision relies on image processing.

Image understanding (feature extraction) may be equally used in both - pure image processing and computer vision.


State Space Observer Control I am attempting to implement observer based state feedback in C and can't figure it out.

Here is pseudocode of the algorithm:

//Initialize Stuff

loop:

//read process variable
y = read(ch1)

//compute control variable
u = -K * xhat

//set process input
write(u)

//update state estimate
xhat = xhat + Ts * (A*xhat + B*u + L*(y-C*xhat))

// repeat loop

xhat is the estimated state
K is the feedback gain
L is the observer gain
y is the system output
u is the system input
Ts is the sampling period
A, B, & C are the state space matrices

Is this how observer feedback is supposed to be implemented?

Here is the actual code I have for the loop:

// Get Process Output from Sensors:
y[0] = sensor1();
y[1] = sensor2();

// Compute Control Variable:
u = 0;
for (i=0; i<4; i++)
    u += K[i]*xhat[i];

// Set Process Input:
write_to_DAC(u);

// Calculate: A*xhat + B*u
for (i=0; i<4; i++) {
    Ax_Bu[i] = 0;
    for (j=0; j<4; j++)
        Ax_Bu[i] += A[i*4+j]*xhat[j];
    Ax_Bu[i] += B[i]*u;
}

// Calculate: y-y_hat
for (i=0; i<2; i++) {
    yhat[i]   = 0;
    y_yhat[i] = 0;
    for (j=0; j<4; j++) {
        yhat[i] +=  C[i*4+j]*xhat[j]; // Calculate yhat = C*xhat
        y_yhat[i] = y[i]-yhat[i];     // Calculate the 'Innovation'
    }
}

// Calculate: L*(y-y_hat)
for (i=0; i<4; i++) {
    L_y_yhat[i] = 0;
    for (j=0; j<2; j++)
        L_y_yhat[i] += L[i*2+j]*y_yhat[j];
}

// Update State Estimate:
for (i=0; i<4; i++)
    xhat[i] += Ts*(Ax_Bu[i] + L_y_yhat[i]);

Assuming $y$ and $u$ are scalars, then the first equation:

$$
u = -K \hat{x}
$$

is just a vector-vector multiplication and can be implemented in C as:

u = 0.0;
for (i=0; i<N; i++)
{
    u += xhat[i]*K[i];
}

assuming that $K$ is $1 \times N$ and that $\hat{x}$ is $N \times 1$.

The multiplication C*xhat can be implemented similarly (because I suspect that $C$ is really $1 \times N$ rather than $2 \times N$).

The multiplications L*(y-C*xhat) and B*u are vector-scalar multiplications which, for the latter, would look like:

Bu = 0.0;
for (i=0; i<N; i++)
{
    Bu += B[i]*u;
}

The multiplication  A*xhat  will be slightly different because $A$ is $N \times N$.


why the ultasound image is noisy by speckle noise? My first question :

I want to know why an image is noisy ?.
why the ultrasound, SAR and Laser... images are noisy by a multiplicative noise not an additive noise?

Another question :

So the dominant noise in ultrasound image is  the Speckle noise, but the confusion in my mind is :

The speckle noise is of multiplicative nature, the mathematical writing is
$$f(x,y)=u(x,u)\eta(x,y)\tag{1}$$
where $f$ is the observed image, $\eta$ the noise and $u$ the dispeckle (recovered) image.

But in some paper (as this) the perturbed images with the speckle noise can write as :
$$f(x,y)=u(x,y)+\sqrt{u(x,y)}\eta(x,y)\tag{2}$$
and $\eta$ is with zero-mean Gaussian variable, the first constraint it must $u>0$ and then the model (2) is characteris of the Gaussian distribution, So this is the case of additive noise

There are plenty of sources of noise in an imaging system. In the domain of sonar systems (which I work in) there is a big problem of reverberation, which is a fancy way of saying "sound from something that we aren't interested in. When using a sonar at sea, you get reverb from the sea surface, sea bed and the water volume. Taking volume reverb as a simple example, this is the water itself (or particles in it) returning sound to the receiver. The best way to relate to this is like driving your car on a foggy night. The sea bed and surface can absorb or reflect sound energy, depending on their composition (e.g. sea state). It's worth noting that this will bugger up your sonar image.

Another overlooked source of noise is in the signal processing itself. Your sensor is not going to be perfect and will quite happily distort both the amplitude and phase of the phenomena its trying to measure. In acoustics, you'll find microphones and hydrophones will be "tuned" for a certain bandwidth, and are quite useless outside of this frequency range. Typically your analogue to digit convertor will have some form of quantisation error with the values its spits out. It simply doesn't have the infinite precision required to represent measurements of the real world. Before we move off the subject, we should be aware that any electronics experiences interference from its  environment (e.g. induction, radiation, poor design). Add all this together and you have an estimate of what's called your "noise floor" and as a designer, you want to both minimise this by design and use other techniques to boost the signal.

Once in the digital world, your sample values are pretty much safe from additional noise. But wait there's more. In acoustics you'll want to electrically steer your array and spatially filter the incoming signal (look in different directions). This is also known as beam forming. The basic principle of beam forming is that you apply delays to the signals you receive and then sum them all together. Using the awesome power of constructive and destructive  interference, incoherent signals cancel out, and coherent signals emphasised. This isn't perfect either because the degree of spatial attenuation varies depending on a number of factors (e.g. array length and element spacing). It's worth pointing out that this aspect of processing features quite a bit of arithmetic on data types that have limited precision. Usually it's not noticeable but its there, trust me.


Poles and zeros map, High pass filter Given the following poles and zeros map I have to identify qualitatively, which filter type is represented (Low-Pass, High-Pass, ...). In the solution it says, that this is a high pass filter.

However, I plotted the amplitude response (see below) and to me this looks more like a notch filter or a band stop.

Where is my mistake in reasoning?

Thank You!

Your pole-zero plot depicts a linear phase, real valued, three coefficient FIR filter, if we assume the zeros are complex conjugates on the unit circle.

The frequency response magnitude is exactly zero at the angle of the zeros on the unit circle. Because of these exact zeros, the log magnitude of the frequency response tends to minus infinity at those frequencies, which gives an impression of notch, however the frequency response is not narrow enough for a typical notch filter. Adding two poles close to those zeros would make it a strong notch filter.

Because of the location of zeros, the frequency response magnitude is larger towards frequencies $\pm \pi$ and smaller at low frequencies and for this reason it is classified as a highpass filter. However it's not a very successful highpass filter either.


Numericaly computing an integral in MATLAB I have the following integral

$$\int_{-\infty}^{\infty} g(t-\tau) \, \mathrm d\tau$$

where $g(\cdot)$ is a known function, e.g., $g(t) = \sin(t)$, and I want to evaluate its value with MATLAB. How should I do this in discrete form? Thanks.

Use integral function.

q = integral(g,xmin,xmax)

For example $g(t) = e^{-t}$

$$q = \int_0^{\infty}e^{-t}\mathrm{d}t = 1$$

g = @(t) exp(-t);
q = integral(g,0,Inf)

More information
Matlab numerical integration


Digital image moments in plain English I'm studying OpenCV and, in computer vision and image processing, people speak of blobs, contours, connected regions, and I sometimes hear the phrase "image moments".

I know of an article on Wikipedia about it, but I think it is too technical. I don't really want to go deep into the math background but I want to know what I'm talking about.

Could someone explain to me what image moments are in plain English?

An image moment is simply a number that characterizes the image, construed as the realization of a spatial random variable. If you've taken any probability class, you should remember the concepts of mean and variance, which are derived from the first and second moments of the random variable (the n'th moment of an r.v. is the expectation of its n'th power). Furthermore, the moments of a random variable collectively prescribe its distribution. In other words, you can reduce a probability distribution to a sequence of numbers, and this is useful when you want to compare distributions numerically.

In image processing, if you want to compare images, you might not want the comparison to be sensitive to minor things like rotation, translation, and scale (since the image remains fundamentally the same). Thus the motivation behind invariant moments you see in the Wikipedia article you cited.


What would be the best solution to synchronize data from two inputs. (Camera and Acceleration sensors) For a video stabilization project, a USB webcam (25 fps) and USB acceleration sensors are being used.

the sensors gives a speed of approximately 50 samples/second.

the derivative of the acceleration aka. 'Jerk' is done in order to localize the start and the end of the shake in the video.

the purpose is to compensate the camera movement in real-time, if it were you how would you :

Match the frame (from camera) with the data of the sensor.
Calculate the number of pixels to compensate from acceleration sensors data.

Thank you.

This is a common problem that most Multimedia systems faces and solves.

Basically, closest to the source, there should be one entity (call it originator) that must collect both the data simultaneously and then stamps the each observation i.e. frames of video and sample of sensor, a time stamp from the local (best suitable clock). After this, originator distributes data in whichever form required.

Essentially a common reference clock for both sources as seen by originator. Any encoding, processing, and consumption of this data must preserve and interpret these time stamps.

This synchronization is of course only as good the source clock itself; there can be jitter and skew in such a clock used by originator. Either the originator follow clock of any primary source (data or video whichever reliable) and stamp it on the other, or derive a local clock using either reliable crystal or GPS or NTP etc.


3 Band Crossover with Linkwitz Riley Filter I would like to write a 3 Band "EQ" with BiQuad Filters but for some reason, I get some phase distortion.

The problem already got discussed here.

I found out that I need to work with Allpass filters in order to compensate the phase response of my filters. Btw, the filters are Linkwitz-Riley 2nd order filters.

But when I try to copy the idea to my project, I still have some distortion.

// low pass calcuation for linkwitz riley filters

float lowCutoff; // the lower cutoff frequency
float highCutoff; // the higher cutoff frequency

double omega = M_PI * cutoffFrequency;
double theta = omega / sampleRate;
double kappa = omega / tan(theta);
double delta = kappa * kappa + omega * omega + 2.0 * kappa * omega;

double a0 = (omega * omega) / delta;
double a1 = 2.0 * a0;
double a2 = a0;

double b1 = (-2.0 * kappa * kappa + 2.0 * omega * omega) / delta;
double b2 = (-2.0 * kappa * omega + kappa * kappa + omega * omega) / delta;

// high pass, only the calculation of a0 to b2 changes
// the intermediate values stay the same

double a0 = (kappa * kappa) / delta;
double a1 = -2.0 * a0;
double a2 = a0;

double b1 = (-2.0 * kappa * kappa + 2.0 * omega * omega) / delta;
double b2 = (-2.0 * kappa * omega + kappa * kappa + omega * omega) / delta;

// allpass filter to compensate phase response:

double q = 0.5;
double alpha = (tan(M_PI*q/sampleRate) - 1) / (tan(M_PI*0.5/sampleRate) + 1);
double beta =   -cos(lowCutoff  * M_PI / sampleRate);

double a0 = -alpha;
double a1 = beta*(1.0 - alpha);
double a2 = 1.0;
double b1 = beta*(1.0 - alpha);;
double b2 = -alpha;;

Since I use 2nd order filters, I also have to use a 2nd order all pass filter, right?
I've read that for those filters and with this calculation, I hav to use q. So q should be the same as the q in the Linkwitz-Riley filter, which should be 0.5 right? (I've read this somewhere, not sure about it.)

Then, my actual processing looks something like this:

// filter names:
// loAPF, lowHPF, ... low: lower cutoffhigh: higher cutoff
// APF: allpass, HPF: highpass, LPF: lowpass

highSample =  lowAPF->tick(highHPF->tick(sample));
midSample =   lowHPF->tick(highLPF->tick(sample));
lowSample =   lowLPF->tick(highLPF->tick(sample));

Does somebody have any idea, why I might have problems with the phase response?

Thanks in advance!

A 2nd-order Linkwitz-Riley crossover needs 1st-order allpass phase compensation.

The Lowpass half of the LR crossover is:

$$\dfrac{1}{s^2+2s+1} $$

The Highpass half of the LR crossover is:

$$\dfrac{s^2}{s^2+2s+1} $$

With even-order LR crossovers, in order to sum-to-allpass one of the two sections has to have inverted polarity. Thus:

$$\dfrac{s^2}{s^2 + 2s +1} - \dfrac{1}{s^2 + 2s +1}= \dfrac{(s-1)(s+1)}{(s+1)(s+1)} = \dfrac{(s-1)}{(s+1)}$$


What is the opposite of sampling called? We call the conversion from a continuous time signal $f(t)$ to a discrete time signal $f_s[k]$ "sampling". Is there a name for the reverse operation, i.e. creating a continuous time signal from a discrete time signal?


Apply Kalman filter to remove process noise, given zero measurement error I have a recording, corrupted by environment noise. The measurement noise is zero. How can I use Kalman filter to remove environment noise? I tried this Matlab code here but if the measurement error covariance is set to zero, the Kalman values are always the same as the measurements. The process noise is not removed at all.

You're misunderstanding how the Kalman filter framework works. Your measurements are a reflection of what the underlying state of the system is at each time instant. If you have zero measurement covariance, then your measurements perfectly reflect the system state, so there's no need for any filter whatsoever. It's not likely that this is truly the case.

Process noise is a representation of uncertainty in the model for how your system's state changes as time goes on. It's not something to be "removed" by the Kalman filter. Instead, it gives you the ability to accommodate higher-order effects that you aren't modeling in your state transition matrix.


Symmetry and periodicity of frequency-shifted discrete Fourier transform The unitary discrete Fourier transform (DFT) of a sequence of numbers $x_n$ to $X_k,$ with integer $0 \le n < N$ and $0 \le k < N,$ can be defined as:
$$X_k = \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x_n e^{-2\pi ikn/N}\tag{1}$$
and the inverse discrete Fourier transform (IDFT) as:
$$x_n = \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X_k e^{2\pi ikn/N}\tag{2}$$
If $x_n$ is modulated (multiplied) by a unit-magnitude zero-phase complex sinusoid $e^{-2\pi ibn/N}$ before DFT, and the IDFT output is demodulated (divided) by the same, then we get another transform pair from the family of generalized discrete Fourier transforms, parameterized by the constant $b:$
$$
\begin{aligned}
X_k(b) &= \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x_ne^{-2\pi ikn/N}e^{-2\pi ibn/N}\\
&= \frac{1}{\sqrt{N}} \sum_{n=0}^{N-1} x_ne^{-2\pi i(k+b)n/N}
\end{aligned}
\tag{3}
$$
$$
\begin{aligned}
x_n &= e^{2\pi ibn/N}\frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X_k(b)\cdot e^{2\pi ikn/N} \\
&= \frac{1}{\sqrt{N}} \sum_{k=0}^{N-1} X_k(b)\cdot e^{2\pi i(k+b)n/N}
\end{aligned}
\tag{4}
$$
Whereas DFT samples frequencies $2\pi k/N,$ the frequency-shifted transform samples frequencies $2\pi (k+b)/N.$ This can be visualized on the Z-plane:
Figure 1. Z-plane representation, showing the unit circle, of the frequencies sampled by A) DFT and B) the frequency-shifted DFT, with $N = 12$ and $b = 1/3.$ The shift $2\pi b/N$ appears as the angle between the real axis and the vector from origin to one of the sampled frequencies.
Question: What are the time-domain symmetry and periodicity properties of such frequency-shifted Fourier transforms when extending $n$ beyond $0\le n<N$ in the inverse transform, and how does this depend on the parameter $b?$
For DFT (or with $b = 0$) the time-domain extension has period $N$ with no time-domain symmetry imposed by the transform.

One version that's useful is $b=\frac{1}{2}$

If you have an input sequence of N real numbers, the output sequence can also be represented as N real numbers: Two for DC and Nyquist and 2*(N/2-1) for N/2-1 complex values. Mixing real and complex numbers makes real time processing a bit more awkward.

Using the transform you suggest with $b=\frac{1}{2}$ results in exactly N/2 complex numbers which allows for more efficient code on SIMD real time processors. $b=\frac{1}{2}$  maintains complex conjugate symmetry. Instead of, say, sampling -10,0,10,20,30, .. you sample at -150,50,50,150,250,... avoiding DC and Nyquist.


Calculating Data Rate using Bandwidth, Transmission Power, Noise Power Spectrum Density and Channel Gain I'm a Software Engineering student researching on Load Balancing in IoT and I need to calculate data rate. However, I don't know how to use these values to attain data rate despite searching on the Internet a lot. Actually, I am not sure whether given values are sufficient or not.

I know that:

Bandwidth (W) = 10 MHz
Noise power density (N0) = -174 dBm/Hz
Path loss model = 128.1 + 37.6log10d, where d refers to the distance in
kilometers.
Transmission power (P) = 23dbm

and the data rate:

$$r =Wlog_2(1 + \frac {P_{tr} H}{N_0W})$$
My questions are:

How should I calculate channel gain (H)? Is it related to path loss
model?
Should channel gain change to the power 2?
Do parameters need any unit conversion?

Let me answer just in terms of the units (the other side-questions have been discussed in the comments already):

Shannons formula is to be interpreted as $W \log_2(1+{\rm SNR})$, where SNR is the signal to noise ratio.
A signal to noise ratio is a ratio of powers: received signal power divided by total noise power at the receiver.
Your received signal power is the transmitted signal power $P$ reduced by what you lose on the way ("path loss").
"Reduced by" means fraction in linear scale (1 Watt transmitted and a path loss of 1000 means you receive 1 Milliwatt).
We often prefer to use logarithmic scales instead since division becomes subtraction in log scale, which is easier. To do this, we define logarithmic measures of power: The logarithmic equivalent of a power in W is $10 \log_{10} \frac{P}{1W}$ and measured in dBW (dB Watt). It means 1 W is the same as 0 dBW and every doubling adds 3 dBs (2 W = 3 dBW, 4 W = 6 dBW, 8 W = 9 dBW and so on). Alternatively, we can use Milliwatt as a reference and get $10\log_{10} \frac{P}{1{\rm mW}}$ in dBm (dB Milliwatt).
You are transmitting 23 dBm and losing $128.1+37.6 \log_{10} \frac{d}{1{\rm km}}$ on the way. Hence, as an example, in 100m distance, you have lost 90.5dB. This means your received power is 23 dBm - 90.5 dB = -67.5 dBm. In linear scale this is $10^{-67.5/10} = 0.177 µW$.
Your thermal noise power density is given as -174 dBm/Hz. This means $10^{-17.4} $W/Hz which at a bandwidth of 10 MHz becomes $10^{-10.4} = 0.04 nW$.
Hence your SNR is approximately $0.177/0.00004 = 4400$. With perfect link adaptation this would give you $\log_2(1+4400) = 12.1$ bits per channel use (That's a crazy high number! I can tell you that you should not expect more than 2-4 bit per channel use in practice.)
At a bandwidth of 10 MHz you could use symbol rates up to $10^7$ channel uses per second.

Just to be safe, a few words of warning when converting this to a data rate. This is risky since a lot of things are missing:

Fading would lead to variations in the SNR which means you'd need to go lower in your modulation order to account for these.
At 10 MHz you'd encounter frequency-selective fading, which requires suitable mechanisms to combat it, e.g., OFDM or SC-FDMA. For OFDM, you then have a frequency-dependent SNR and need to apply bit loading to make proper use of all the subcarriers. IOT devices would do something simpler like SC-FDMA but probably not at 10 MHz.
Reliable communication will require coding, which comes with some overhead.
You cannot fully use the 10 MHz, due to spectral leakage of whatever pulse shape you are using, you'd need to plan for some guard bands.
The noise calculation is overly optimistic since no noise figure was included. Typically, you would need to increase your receiver noise power by the noise figure of your LNA. IOT devices may use simple LNAs with higher noise figures.
In IOT you have multiple access problems to solve, i.e., the devices need to share the spectrum. This comes at some loss, since there will always be times where the spectrum is underutilized and you need to leave room for the initial access of new devices.

And that is still an incomplete list :)


What is the difference between Lumped and Distributed systems? What are the salient differences between Lumped and Distributed systems? In what contexts are distributed systems the appropriate model and in what are lumped systems the appropriate model?

Also, Lumped systems are said to be described by ordinary differential equations while the latter is said to be described by partial differential equations. Can someone explain why?

The elements building a lumped system are thought of being concentrated at singular points in space. The classical example is an electrical circuit with passive elements like resistor, inductance and capacitor. The physical quantities current and voltage are functions of time (only). E. g. the current at a capacitor with capacity $C$ is given by
$$
i(t) = C\frac{\mathrm d v(t)}{\mathrm d t}
$$
Where $C$ is a constant (and so are $R$ and $L$). This leads to ordinary differential equations.

In contrast, the elements in distributed systems are thought of being distributed in space, so that physical quantities depend on both time and space. The classical example is the electrical line where inductance, capacity and resistance are not constant but functions of length $x$. This leads to partial derivatives of $i(t,x)$ and $v(t,x)$ in $t$ and $x$.


integration property of fourier series Please help me sort this issue out.
The integration property in Fourier series is as follows:

So, for proving the above property, i took this approach:

This is where my doubt is. Some books and websites just put upper limit (ignoring the lower limit) and compare with (1) to conclude that

However, when the lower integral limit is substituted, we get

This value cannot be evaluated. Though some websites say that this corresponds to the "initial value"(?) and that is zero... etc., how does one justify this mathematically?
Any help is appreciated. Thank you.

Note that the antiderivative of a function is only defined up to a constant. Furthermore, note that if you integrate a periodic function, the result is not necessarily periodic. Let

$$x(t)=\sum_{k=-\infty}^{\infty}a_ke^{jk\omega_0t}\tag{1}$$

Now we integrate $(1)$ with an arbitrary lower integration limit $t_0$. I'll say more about that later.

$$\begin{align}y(t)&=\int_{t_0}^tx(\tau)d\tau\\&=\sum_{k=-\infty}^{\infty}a_k\int_{t_0}^te^{jk\omega_0\tau}d\tau\\&=a_0(t-t_0)+\sum_{k=-\infty\\k\neq 0}^{\infty}\frac{a_k}{jk\omega_0}e^{jk\omega_0t}-\sum_{k=-\infty\\k\neq 0}^{\infty}\frac{a_k}{jk\omega_0}e^{jk\omega_0t_0}\tag{2}\end{align}$$

Clearly, if $a_0\neq 0$, $y(t)$ isn't periodic, so we can't hope to obtain a formula for its Fourier coefficients. So in order to end up with a periodic function $y(t)$, we require that $a_0=0$.

If $a_0=0$, we get from $(2)$

$$y(t)=\sum_{k=-\infty}^{\infty}b_ke^{jk\omega_0t}\tag{3}$$

with

$$b_k=\begin{cases}\displaystyle\frac{a_k}{jk\omega_0},&k\neq 0\\\displaystyle-\sum_{l=-\infty\\l\neq 0}^{\infty}\frac{a_l}{jl\omega_0}e^{jl\omega_0t_0},&k=0\end{cases}\tag{4}$$

So in general the antiderivative $y(t)$ has a non-zero DC component, which depends on the choice of the lower integration limit $t_0$.

The antiderivative

$$\tilde{y}(t)=\sum_{k=-\infty\\k\neq 0}^{\infty}\frac{a_k}{jk\omega_0}e^{jk\omega_0t}\tag{5}$$

is the specific antiderivative of $x(t)$ that has a zero DC Fourier coefficient, and this is the one that is meant by the formula that is usually given in textbooks.

Consequently, the expression

$$y(t)=\int_{-\infty}^tx(\tau)d\tau\tag{6}$$

for the antiderivative of $x(t)$ with Fourier coefficients $b_k=a_k/(jk\omega_0)$, $k\neq 0$, and $b_0=0$,
is at least inaccurate; it is in fact a sloppy way of expressing the specific antiderivative $\tilde{y}(t)$ given by $(5)$, i.e., the one with a zero DC coefficient.


Short length windowed sinc filter I have a 50-60Hz signal being sampled at a frequency of 25 kHz & I'm trying to filter all frequencies beyond 70 Hz with a reasonable steep transition band (approx 4 Hz) & good attenuation (~40+ dB).

For this, I tried a Windowed sinc filter using the following parameters :-

Cutoff frequency, fc = 70 / 25000   = 0.00280

Bandwidth,        BW = 4  / 25000   = 0.00016

Filter length,    M  = 4  / 0.00016 = 25000

The problem here is the filter length. The value of 25000 means a window that's 1 second wide while my target is typically around 500-1000 samples, i.e. 1-2 cycles at 50 Hz.

I can't increase the bandwidth & I'm not sure what other options I have to meet this filtering requirement. So, my questions are,

Is the windowed sinc not a suitable filter for this application?
Considering that the sampling frequency is so high, shouldn't it be quite easy to achieve this type of filtering with a small filter kernel?
Is there any other filter that I can use for this? My requirement is to offer good attenuation to all signals beyond 70 Hz using a small window. The choice of the filter isn't of much concern.

The following is all assuming we are interested in Finite Impulse Response (FIR) filter design.

The cutoff frequency is not really the determining factor in the filter length. The 3 main things that determine the filter length are: the amount of passband ripple, the amount of stopband attenuation and the relative size of the transition band, i.e., the size of the band going from the passband to the stopband. Cutoff frequency will affect the design if it is really close to 0 Hz or half your sampling frequency.

The smaller the size of the transition band the longer your filter is going to be. To specifically answer your questions.

The sinc filter design is largely by trial an error to achieve your desired specifications. You play with the length of filter, the type of window used, and the filter cut-off to try to meet your specifications. So it may not be the most suitable in terms of trying to meet your spec. That said, it will give you an idea of the size of filter that would be required. Other design techniques are likely to require filters of similar lengths.

Regarding the transition width and sampling frequency, as I mentioned above it is actually the opposite. The smaller the transition width, the longer your filter is going to need to be.

One alternative, is to actually downsample your signal to reduce your sampling frequency. Since you are attenuating all signals above 70 Hz, why not just use decimation? The advantages of this technique are:

Reduced computational load
The filters used can be much smaller because they can have larger transition widths.
If really needed, you can upsample to get the data back to the original sampling rate using the same filters you used for the downsampling.

The main disadvantage is that your overall design is a bit more complicated. A good reference for this is "Multirate Signal Processing" by Crochiere and Rabiner.

If you are looking for the smallest filter lengths, then using the Parks–McClellan algorithm based on the Remez exchange technique has the shortest filter lengths for meeting passband ripple and stopband attenuation levels.


Gesture recognition from analog sensor value Not sure if this is the right place for my question but let me give it a try.
I have a very simple analog sensor (piezo transducer) which is attached to a microcontroller. I'm using this sensor to detect gestures (hand movements).
The mCU sends the sensor data to my PC where I'm running a simple application that plots the data to a graph. In this graph the human eye and brain can recognize easily which gesture is present on the graph, but, of course, the computer has to recognize it as well.

I'm looking for an algorithm or at least a starting point I can detect these gestures. I'm very new to this field so any advice or information would be a great help.

Here you can see a sample screenshot of the graph with two different gestures:

Feel free to set the correct tags for this question.

You can use a matched filter. Average all the gesture 1's and that's your filter. Do the same for however many gestures you have and you will effectively have a bunch of filters. Run the data through those filters in parallel and then you will need a decision threshold that indicates yes/no if that gesture is detected.

It should look something like this:


Are there any signal processing libraries that implement TD-PSOLA? For such a popular and old algorithm, I've had surprisingly little luck finding implementations. I've got it working in MATLAB, but am curious if any other languages have libraries containing it (python?). I've googled extensively to no avail.

When I'm looking for source code, I generally search the open repositories rather than Google.

Try this one in C++.

However, there don't seem to be many on github, and sourceforge has nothing of interest.


What is the most efficient way to perform Gaussian blurs for varying sizes of Gaussian kernels in an image? I have detected corner points of images using Harris corner detector in GLSL and now, I intend on using Fast Retina Keypoint(FREAK) as keypoint descriptor to later match the keypoints between images.

FREAK uses circular sampling pattern, where each circles represent Gaussian kernels applied to the corresponding sampling points.

And, the size of Gaussian kernel increase exponentially for points farther away from central keypoint. I previously used nine hit Gaussian blur as mentioned here. So, what is the most efficient way to perform Gaussian blurs for a number of Gaussian kernels with varying sizes, which also overlap themselves, in an image?

I'm not really used to thinking in texture fetches as I usually code high level non-hardware-specific signal processing, but low pass filtering also removes information you can couple it with sub-sampling which can give further gains. It is not clear to me if you have already done this.

In practice this means that after the first low-pass filtering, you can disregard $k$ out of $N$ of the samples in the subsequent low-pass filtering.

Simple example is the $2^n \times 2^n$ averaging-filter, first we can do 2x2 mean values and store. Then for the next pass, we only need to use the output of one of each $2 \times 2$ blocks. So we only need consider $1$ in $4^k$ data points at step $k$.

A 8x8 block this would mean:

$$I_1 = \left[\begin{array}{cccccccc}
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1\\
1&1&1&1&1&1&1&1
\end{array}\right]$$
First iteration, all are involved in 2x2 filtering. But they are only stored at the pixels which are ones here (so we only need to calculate 1/4 of them):
$$I_2 = \left[\begin{array}{cccccccc}
1&0&1&0&1&0&1&0\\
0&0&0&0&0&0&0&0\\
1&0&1&0&1&0&1&0\\
0&0&0&0&0&0&0&0\\
1&0&1&0&1&0&1&0\\
0&0&0&0&0&0&0&0\\
1&0&1&0&1&0&1&0\\
0&0&0&0&0&0&0&0
\end{array}\right]$$
Next step only involves the ones in the previous one, now our image is 4x4 instead of 8x8, so calculations are reduced a factor of 4. But we only need to store the results at pixels:
$$I_3 = \left[\begin{array}{cccccccc}
1&0&0&0&1&0&0&0\\
0&0&0&0&0&0&0&0\\
0&0&0&0&0&0&0&0\\
0&0&0&0&0&0&0&0\\
1&0&0&0&1&0&0&0\\
0&0&0&0&0&0&0&0\\
0&0&0&0&0&0&0&0\\
0&0&0&0&0&0&0&0
\end{array}\right]$$
We will end up with only one out-value, in this case upper left. And here's where the interpolations come into play. We can interpolate between these upper left results of neighboring "blocks". The interpolation kernel can be dependent on the low-pass filter, but I assume that this is where the hardware shines with it's strengths, having really fast ways to do this for us?

If you are interested in theory filter-banks, sub-band coding and maybe wavelets are interesting to read up on. The example I brought up here is the implementation low-pass filter of the basic but popular discrete Haar-wavelet. Gaussian pyramids in computer vision may be even more relevant. There's much to read about it so it's difficult to cram it all into just one answer.


In framing of audio samples', what is need of frame shift while giving frame size?? In framing of audio samples' in audio feature extraction, what is need of frame shift while giving frame size??
i.e. frame size  = 20ms,
     frame shift = 10ms.
Rather then shifting/overlapping why can't we use continuous frames then of overlapped??

By performing the windowing with overlap we are artificially increasing our time resolution (larger granularity of features in time). This is especially useful when frame duration is long (bad time resolution, very good frequency resolution), thus yielding kind of extra 'time resolution'.
Usually no one is using the rectangular window, but other types such as Hamming, Hann, etc. As you probably know, windows tend to zero samples at both of their ends. By doing data you are kind of loosing some data. Let's imagine that you have some very short event you want to detect, at point where two windows are joined. If there is no overlap, then these samples will be almost removed. But if you perform 50% of overlap, there will be another window, right in the centre, and you won't lose anything.