Calculate harmonics using DFT from real points I have a real data of 144 points, when I perform a 144-point DFT on this data, I get $X$ with real and complex values. I want to calculate harmonics using these $X$'s.

The $X[0]$ and $X[72]$, added together and divided by 144, would give me the DC component?
And can I just use the next 71 $X$'s and their conjugates (Euler's identity), to calculate the harmonics? I believe it's $(X[i] + X^*[i])/144 $?

Sorry about that, but I have this problem to solve, I'll try and explain it and what I'm supposed to do with it.
Thanks for being patient.

I have a data set of 144 points. I want to express this in terms of a DC component and it's harmonics. So, I was asked to perform a DFT operation on this and get the harmonics that would represent this data in frequency domain. I performed a 144 point DFT on this data and got 144 $X[k]$'s with real and imaginary parts. As, per my understanding these 144 $X[k]$'s ($X[0]$ being the DC) represent the time domain signal in frequency domain, but I'm still being asked for the DC(apparently ($X[0]+X[72])/2$) and the 71 harmonics, which I'm not sure how to go about.

I apologize if this still doesn't makes sense, but it's what I'm supposed to do, probably you can give me some reference where I can get such concept or anything related.

The DFT will express your time sequence as a weighted some of complex exponential (basically a set of orthogonal functions, if that's a helpful concept for you). DC is simply a special case of a complex exponential with the frequency being 0.

Let's assume a sample rate of 1440 Hz (to make the math simple). Then the DFT coefficients mean the following:

X[0]: amplitude at 0 Hz
X[1]: amplitude & phase at 10 Hz
X[2]: amplitude & phase at 20 Hz
...
X[71]: amplitude and phase at 710 Hz
X[72]: amplitude at 720 Hz (Nyquist frequency)

Nyquist is a little bit of an odd-ball. It's a real number so the phase is 0. The value should also be small or close to 0, otherwise you have potentially aliasing.

Since your input is real, you have complex conjugate symmetry as follows

X[73] = X[71]';
X[74] = X[70]';
....
X[143] = X[1]';

Just a quick comment: (X[0] + X[72)/2 is NOT the DC value. X[0] is.


ADC and Matched Filtering A basic theorem in communications is the matched filter maximizes the SNR at sampling. I'm a little confused on how this relates to discrete time systems and sampling rate.

Normally if you sample at Nyquist rate it's "good enough" but how does this interact with the discrete time correlator? Does sampling at a higher rate benefit the system? Does up-sampling the signal and passing it through the correlator yield the same benefit?

I understand the idea of Noise being White and the output of white noised passed through the correlator limits its PSD but lets say your analog signal is pre-filtered to remove bands outside of some region and your sample rate meets the Nyquist rate for this bandwidth. Is there any benefit to over-sampling and then discrete time matched filtering?

Normally if you sample at Nyquist rate it's "good enough" but how does this interact with the discrete time correlator?

It interacts very well! See the following:

  Does sampling at a higher rate benefit the system?

You mean oversampling.

If you filter your signal to the original signals bandwidth, which is now smaller than the Nyquist rate. You thus get a component that is correlated between samples (which is your signal) and one that isn't (which is your noise). By averaging / downsampling your signal, you increase the variance of your signal component more than that of your noise component – you get an SNR gain.

  Does up-sampling the signal and passing it through the correlator yield the same benefit?

No, because upsampling doesn't "uncorrelate" the noise that you already digitize.

  Is there any benefit to over-sampling and then discrete time matched filtering?

Yes, as described above, the oversampling allows you to increase your SNR. Think of it this way:

Assume that if critically sampled, your pulse length is $N$; thus, your MF does a correlation over $N$ samples.

You oversample your analog signal by a factor of $f>1$. You interpolate your MF by the same factor, so that the pulse and the filter are now $Nf$ in length.

While the power of the samples stays the same, the amplitude of the signal part in the  (digital) correlator output has been scaled by a factor of $f$, which means the power  has been scaled by $f^2$.

For the noise component in the correlator output: White Noise samples are uncorrelated, and we know that thus the variance of their sum (which is what a correlator is, weighted sum) is the sum of their variances. Since variance is power for zero-mean noise, that means we get a noise power increase of $f$.

Thus, oversampling by $f$ gets us an SNR increase of $\frac{f^2}f$ :)


Radix2 fft - zero padding output size With zeros padding the FFT radix2 has different input size than output. How to deal with that? I see two solutions but not sure which one is better/appropriate. Maybe you know some other and better solutions?

Let me explain:

let’s sey I have signal buffer size 400, and I want to perform radix2 fft with zeros padding.
So the nearest $2^L$ is 512. So I need to add 112 zeros for each buffer.
Then in the fact the output buffer is 512. But I need 400, the same as input.

So first solution I see is to provide some FIFO algirithm for output signal.

And the second solution I see is to shrink my output from 512 to 400, and spread all 512 bins evenly on 400 bins. But then adjacent bins will partly overlap each other. So I need to avarage them in some way. HOW? And is it at all solution that would give me nice results, for example for pitch shifting?

For any help thanks in advance.

Note the convention of notation: "$X[\cdot]$" is not the same as "$X(\cdot)$" (but they are defined solely in terms of the same "$x[\cdot]$").  The square brackets [] allow only discrete integer arguments, the round parenths () allow continuous real or complex arguments.

Discrete-Time Fourier Transform (DTFT):

$$ \mathscr{F}\big\{x[n]\big\} \triangleq X\big( e^{j\omega} \big) = \sum\limits_{n=-\infty}^{\infty} x[n] \, e^{-j\omega n} $$

This is literally the double-sided $\mathcal{Z}$-transform of $x[n]$ evaluated at $z=e^{j \omega}$.

If your original $x[n]$ has non-zero length of $N$, then DTFT is:

$$ X\big( e^{j\omega} \big) = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j\omega n} $$

The length-$N$ Discrete Fourier Transform (DFT) of $x[n]$ is

$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \, e^{-j 2 \pi n k /N} $$

The DTFT is a mathematical operation for the blackboard (or whiteboard or on paper).  The DFT is what computers do when they perform the Fast Fourier Transform, (FFT).

And the DFT is the DTFT of this finite-length sequence evaluated at $N$ equally-spaced frequencies between $0$ and $2\pi$.

$$ X[k] = X\big( e^{j\omega} \big) \bigg|_{\omega = 2 \pi k/N} $$

So now, what if you zero-pad $x[n]$ to a longer length $\tilde{N}>N$?

$$ \tilde{x}[n] = \begin{cases}
x[n] \qquad & 0 \le n < N \\
\\
0 \qquad & N \le n < \tilde{N} \\
\end{cases} $$

Then the DTFT of this is:

$$ \mathscr{F}\big\{\tilde{x}[n]\big\} = \tilde{X}\big( e^{j\omega} \big) = \sum\limits_{n=-\infty}^{\infty} \tilde{x}[n] \, e^{-j\omega n} $$

But the funny thing is that the two DTFTs are the same

$$ \tilde{X}\big( e^{j\omega} \big) = X\big( e^{j\omega} \big) $$

because there is no difference between $x[n]$ and $\tilde{x}[n]$ in the zero-extended DTFT.

But the DFTs come out a little different (only because $\tilde{N} \ne N$).

$$\begin{align}
\tilde{X}[k] &= \sum\limits_{n=0}^{\tilde{N}-1} \tilde{x}[n] \, e^{-j 2 \pi n k /\tilde{N}} \\
\\
&= \sum\limits_{n=0}^{N-1} \tilde{x}[n] \, e^{-j 2 \pi n k /\tilde{N}} \\
\end{align}$$

(Note the change of $\tilde{N}$ to $N$ in the second summation is because $\tilde{x}[n]=0$ for $N \le n < \tilde{N}$.  That has an interesting consequence below.)

Also this:

$$ \tilde{X}[k] = \tilde{X}\big( e^{j\omega} \big) \bigg|_{\omega = 2 \pi k/\tilde{N}} $$

Now it turns out that the inverse DFT is similar for $\tilde{x}[n]$

$$ \tilde{x}[n] = \tfrac1{\tilde{N}} \sum\limits_{k=0}^{\tilde{N}-1} \tilde{X}[k] \, e^{+j 2 \pi n k /\tilde{N}} $$

and the DTFT can be evaluated in terms of this expression for $\tilde{x}[n]$.

$$\begin{align}
 \tilde{X}\big( e^{j\omega} \big) &= \sum\limits_{n=0}^{\tilde{N}-1} \tilde{x}[n] \, e^{-j\omega n} \\
\\
&= \sum\limits_{n=0}^{N-1} \tilde{x}[n] \, e^{-j\omega n} \\
\\
 &= \sum\limits_{n=0}^{N-1} \tfrac1{\tilde{N}} \sum\limits_{ k=0}^{\tilde{N}-1} \tilde{X}[k] \, e^{+j 2 \pi n k /\tilde{N}} \, e^{-j\omega n} \\
\\
 &= \sum\limits_{k=0}^{\tilde{N}-1} \tilde{X}[k] \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi n k /\tilde{N}} \, e^{-j\omega n} \\
\end{align}$$

The DFT of $x[n]$ you seek is

$$\begin{align}
X[k] &= X\big( e^{j\omega} \big) \bigg|_{\omega = 2 \pi k/N} \\
\\
&= \tilde{X}\big( e^{j\omega} \big) \bigg|_{\omega = 2 \pi k/N} \\
\\
&= \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi n \tilde{k} /\tilde{N}} \, e^{-j\omega n} \bigg|_{\omega = 2 \pi k/N} \\
\\
&= \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi n \tilde{k} /\tilde{N}} \, e^{-j2 \pi n k/N} \\
\\
&= \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)n} \\
\end{align}$$

Note the $\tilde{k}$ is just a different "$k$" as a token counting variable in the summation.

There is this exponential summation identity in closed-form:

$$ \sum\limits_{m=0}^{M-1} a^m = \frac{a^M-1}{a-1}$$

so

$$\begin{align}
\tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)n}  &= \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  \big(e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)}\big)^n \\
\\
&= \frac{\big(e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)}\big)^N - 1}{\tilde{N} \big(e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)} - 1 \big)} \\
\\
&= \frac{e^{+j \pi (\tilde{k}N/\tilde{N} - k)}\big( e^{+j \pi (\tilde{k}N/\tilde{N} - k)} - e^{-j \pi (\tilde{k}N/\tilde{N} - k)} \big)}{\tilde{N} e^{+j \pi (\tilde{k}/\tilde{N} - k/N)}\big(e^{+j \pi (\tilde{k}/\tilde{N} - k/N)} - e^{-j \pi (\tilde{k}/\tilde{N} - k/N)} \big)} \\
\\
&= \frac{e^{+j \pi (\tilde{k}N/\tilde{N} - k)}e^{-j \pi (\tilde{k}/\tilde{N} - k/N)}\big( e^{+j \pi (\tilde{k}N/\tilde{N} - k)} - e^{-j \pi (\tilde{k}N/\tilde{N} - k)} \big)/(2j)}{\tilde{N} \big(e^{+j \pi (\tilde{k}/\tilde{N} - k/N)} - e^{-j \pi (\tilde{k}/\tilde{N} - k/N)} \big)/(2j)} \\
\\
&= \frac{e^{+j \pi (N-1)(\tilde{k}/\tilde{N} - k/N)} \sin\big( \pi N (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}{\tilde{N} \sin\big( \pi (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)} \\
\end{align}$$

This shows the DFT of the original data (at the original size of $N$) can be extracted from the DFT of the zero-padded data:

$$\begin{align}
X[k] &= \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \tfrac1{\tilde{N}} \sum\limits_{n=0}^{N-1}  e^{+j 2 \pi (\tilde{k}/\tilde{N} - k/N)n} \\
\\
&= \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \left( e^{+j \pi (N-1)(\tilde{k}/\tilde{N} - k/N)} \frac{\sin\big( \pi N (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}{\tilde{N} \sin\big( \pi (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)} \right)  \\
\end{align}$$

It ain't purty, but I think this is the simplest I can make it.  This is one form of expressing the Dirichlet kernel.

Edit note:

I made a completely legitimate change of $\tilde{N}$ to $N$ (because $\tilde{x}[n]=0$ for $N \le n < \tilde{N}$) and this results in some changes in the answer.  With $\tilde{N}$ unchanged the previous answer (which is just as correct) was:

$$ X[k] = \sum\limits_{ \tilde{k}=0}^{\tilde{N}-1} \tilde{X}[\tilde{k}] \left( e^{+j \pi (\tilde{N}-1)(\tilde{k}/\tilde{N} - k/N)}\frac{\sin\big( \pi \tilde{N}(\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}{\tilde{N} \sin\big( \pi (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)} \right) $$

They both have to be correct.  But clearly these two expressions are not equal:

$$ e^{+j \pi (\tilde{N}-1)(\tilde{k}/\tilde{N} - k/N)}\frac{\sin\big( \pi \tilde{N}(\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}{\tilde{N} \sin\big( \pi (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}
\ne e^{+j \pi (N-1)(\tilde{k}/\tilde{N} - k/N)}\frac{\sin\big( \pi N (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)}{\tilde{N} \sin\big( \pi (\tfrac{\tilde{k}}{\tilde{N}} - \tfrac{k}{N}) \big)} $$

So how is it that the summations (with the very same $\tilde{X}[\tilde{k}]$ coefficients) add up to be the same?

I would be very interested in anyone's answer to that.


Explanation a flat fading channel in case MIMO OFDM I am trying to understand flat fading and frequency selective channels.

As I understood, in the case of flat fading channels  MIMO channel elements, $h_{i,j}$, are random variables.

Which assumption can I make, if the channel flat fading of frequency selectiv for channel matrix?

Actually, I am reading about MIMO OFDM. In all paper which I have read, frequency selective channel model is given for MIMO OFDM. But in my task I need to assume a flat fading. Can i do such assumption?

Frequency-flat fading essentially means that the channel's transfer function is (approximately) constant in your frequency band of interest, which means you can treat it as if it would not depend on frequency at all. This means two things:

Your MIMO transmission model is simply $\mathbf y(t) = \mathbf H \mathbf x(t) + \mathbf w(t)$ with a transmit signal $\mathbf x(t)$, a MIMO channel matrix $\mathbf H$, a received signal $\mathbf y(t)$ and your additive noise $\mathbf w(t)$. Note that $\mathbf x(t)$ needs to be narrow-band for this to make sense.
Your wouldn't use OFDM in this case. It just doesn't make sense: Every subcarrier would see exactly the same channel and you would gain nothing from the OFDM pre/post processing.

As for modelling the elements of your channel matrix, you can use statistical models that assume some random distribution that the elements are drawn from but be wary of the assumptions those are derived under. For instance, if you assume dense multipath from all directions, you end up with a Rayleigh distribution (resembling a case where you have rich scattering and no line of sight component). For cases with line of sight you might want to consider a Rician distribution. This is fine for not too high frequencies. For mmWave bands, this would be a bad model, since there the multipath is typically quite sparse and your channel matrix tends to be low rank.


Confidence Interval For Frequency of Peak in Averaged Periodogram? I am averaging periodograms from sets of ~15-25 experimentally measured timeseries all of the same length.  To do this I am just summing the power in each frequency bin over all periodograms and then dividing by the total power over all the periodograms (I believe this is similar to Welch's Method but here I have many timeseries instead of a single long timeseries).

In many of these averaged periodograms I see a distinct peak that is called as statistically significant by Fisher's G-test.  What is a good way to get a confidence interval or other interval estimate (e.g. standard-error of the mean) for the frequency of this peak?  Is there some method that takes into account both sample size and width of the peak (maybe width at half-max)?

Add note below to provide more detail:

These are real timeseries data from measurements of gene expression in cells (each timeseries is the gene expression in an individual cell sampled every 6 minutes for ~10 hrs).  I have been averaging the individual periodograms  to create an averaged periodogram (like Bartlett's method) and want to be able to determine two things:

i.) Is a power at a frequency in the averaged periodogram statistically significant (i.e. unlikely to be due to chance)? This is what I am currently using Fisher's G-Test for but if there are better methods I am open to trying those (provided I can understand them and they are relatively easy to use).

ii.)  What is an interval estimate for the frequency of a statistically significant power in the averaged periodogram?  The bin with the greatest power in these averaged periodograms is typically a "peak" (when a line is plotted connecting the power values across bins) with neighboring bins also showing some increased power relative to other parts of the periodogram.

I am averaging periodograms from sets of ~15-25 experimentally measured timeseries all of the same length. To do this I am just summing the power in each frequency bin over all periodograms and then dividing by the total power over all the periodograms...

So, for a given waveform, your time series dataset is $x_i$, where $i$ denotes the $i^{th}$ acquisition, let's call it $i \in [1 \ldots N_{acq}]$. You then derive the periodogram as $X_i = |\mathscr{F(x_i \otimes x_i)}|$. Now, for a given harmonic $k$, each of the $X_i[k], i \in [1 \ldots N_{acq}]$ gives you some point estimate for that $k^{th}$ harmonic.

Following this, what you do is to generate some $\bar{X} = \frac{1}{N_{acq}} \cdot \sum_{i=1}^{N_{acq}}X_i$. So, now, your $\bar{X}$ is the averaged periodogram.

As it seems, the next thing that you want to do is to answer the question "Is this periodogram obtained by chance or is there something in there?"

  In many of these averaged periodograms I see a distinct peak that is called as statistically significant by Fisher's G-test.

OK, so, for this one, you need a test. You are going for the G-Test. The G-Test works over "counts". The value of the frequency bin would be represented as a random variable and we would be asking the question, does this look totally like noise OR is there a significant difference between the values of a frequency bin between experimental conditions (?).

So, the question here is, are your variables so much quantised as to be considered "categorical" and justify the use of the G-Test? (Probably not (?)).

The other point which is not entirely clear is "what" do you want to work on?

At the moment, given a $k^{th}$ harmonic, you work on a dataset of $N_{acq}$ values, which are indexed by $i$. But, what you seem to be wanting to do is work across the $k$s by locating "peaks".

So, the question here is, do you want to be able to say "Over so many experiments, the power of this particular harmonic is some $\bar{X}[k] \pm \sigma$" or do you want to say "Over so many experiments, the frequency is some $f_{critical} \pm \sigma$"?

If you want to work over the $i$s (given some $k$) and therefore make a statement on the validity of the harmonic itself: Assume that frequency bins are independent with each other and go for a "straightforward" t-test. This can either be 1-sample to work out significance or 2-sample to work out significance across two conditions. In using the t-test, you would also be assuming that the distribution of the $k^{th}$ harmonic is distributed normally, in which case you can also apply confidence intervals in the typical way.

If you want to work over the $k$'s, (having averaged in the "$i$ direction") and therefore make a statement on the frequency itself: Then, you cannot assume that the bins are independent with each other (which, they are not) and you would be seeking a model that describes $\bar{X}$ as a whole. One option could be kernel density estimation. Because it would explain $\bar{X}$ like $\hat{X} = \sum_{k=1}^{N_{kernels}} \mathcal{N(\mu_m,\sigma_m^2)}$, it would give you an estimate for the "bump" that is created by the frequency "lobe" at frequency $\mu_m$ with a spread of $\sigma_m^2$. In this case, again, because you have those $\mu_m,\sigma_m^2$ pairs, you could use them to apply confidence intervals.

As a last note, have you tried plotting boxplots of the spectra over the two conditions? Is there overlap in the frequency range you are interested in?

Hope this helps


Image Segmentation - Berkeley Dataset I've been given an image from Berkeley dataset to segmentate. I am an undergrad student and so far I threw everything i know to this thing. But the colors are extremely close. I tried;

Clustering
Edge Detection
Histogram Thresholding

I need suggestions. Which path should I take? I need to segmentate the animals from this picture but I was only able to remove the green leafs.

You can try grabcut algorithm. Check this link:
Interactive Foreground Extraction using GrabCut Algorithm


Signal distortion under non-linear channel The channel is given to have these characteristics, and its required to study the affect on Frequency Modulated signal and its demodulation.
Since there are no sharp corners (like during clipping) will there be a significant difference in the frequencies? Can a regular high-pass filter be used still to demodulate the FM signal.

I think this will help you: Think about how the "frequency" is given by the zero crossings of the waveform, and then ask if the transfer function you show (and for that matter even if it did have sharp corners and clipping), would that have any effect possibly on the zero-crossings? Can we determine the frequency completely and accurately from the zero crossings for a real waveform?
Also consider this: All noise can be decomposed into AM and PM components, meaning Amplitude modulation and Phase Modulation. (Frequency Modulation and Phase Modulation are similar: If you have FM then you have PM since they are related by an integral-- a change in phase over time is Frequency). So if you have a signal that has AM and FM components, and you hard limit it, you remove all the AM components and you would only be left with FM. Here is the thing to think about then: What if you have a FM signal and hard limit it (which is the sharp corners and clipping you mention)-- it would remove AM if it existed, but how would it effect the FM?
For more mathematical thoroughness, explore how ANY pattern that repeats in time over period T can exist in frequency only at frequency n/T for n as the set of all integers. Would the transfer function from x to y that you show possibly change T for any given waveform that repeats over T?


Is there a difference in bicubic interpolation of images vs patches? I need to extract low resolution patches from an image. For example, 32x32 patches from a 500x500 image.

Method One

Down res 500x500 image (by bicubic interpolation, both down res and up res).
Extract 32x32 patches.

Method Two

Extract 32x32 patches.
Down res (by bicubic interpolation, both down res and up res) 32x32 patches.

Will the patches extracted be identical? If they aren't, how different will they be?

The short answer is no. The order of operations doesn't matter and the final result would be similar within reasonable approximation errors.

The long answer is:

A given image $I$ has a pixel resolution of $M \times N$ pixels and a pixel size of $P \times Q$mm. Therefore, the image physical size is $M \cdot P \times N \cdot Q$mm.

If you downsample that (with any interpolation method) by some factor $d \in \left(0 \ldots 1\right)$, then its new physical size will be $d \cdot M \cdot P \times d \cdot N \cdot Q$. Obviously, here $d \cdot M \le 32$ and $d \cdot N \le 32$.

By the way, since at this point we are talking about reduction in physical size, this means that some of the original pixels, some of the original actual information content, had to be thrown away

If you then select the $32 \times 32$ patch, you are selecting some "space" out of the "shrunk" image, some of which information content has now been thrown away. ($32 \cdot d \cdot M \cdot P$)

If you were to select a $32 \times 32$ patch out of the original $M \cdot P \times N \cdot Q$ image, then you select a patch at full resolution.

Once you have selected the patch, the rest of the image is "thrown away". Any downsampling of that patch will produce a new patch whose dimensions (and actual information content) will be lower than $32 \times 32$. If you are thinking at this point "Oh, its alright, I will cut out a $32 \times 32$ patch, downsample it to (for instance) $12 \times 12$ and then upsample that patch back to $32 \times 32$" then, in terms of actual information content you are not gaining anything and you are back to the first case.

In fact, practically, because of the double interpolation (one to downsample and then one to upsample), this second case will have more opportunity for numerical errors but in terms of actual information content, actual spatial resolution, you don't gain anything.

The "dangerous" bit here is the downsampling. Once that takes place, that is where actual information is thrown away. This is why, when you are trying to resample a signal by a non-integer factor $\frac{u}{v}$, you first upsample by $u$ (where no new information content is added but the same information is distributed over larger space) and then downsample by $v$ (where actual information get's limited).

In your case, you have two different routes that both include downsampling by the same factor at some point. Their end result will be similar within reasonable approximations afforded by the interpolation steps.

Hope this helps.


Generating a high SNR sinusoide with lut on a DSP I'm generating a sinus by using a lut method on a DSPic33f. My sample rate is 48 kHz, so I saved 12000 of the first value (unsigned int, 16 bits) and use trigonometric formulas to calculate the other parts of the sinusoid. I can change the frequency of this sine by turning a rotary encoder. Every 2 ms, a timer interrupt read the value of the rotary encoder, and the value (frequency in Hz) will be the new index to read my lut.

A bit of the code for you to understand :

ANGLE += FREQUENCY_VALUE;
if(ANGLE < SINE_TABLE_LEN)
{
    index  = ANGLE;
    sign = 1;
}
else if((ANGLE >= SINE_TABLE_LEN) && (ANGLE < SINE_TABLE_LEN2))
{
    index = SINE_TABLE_LEN2 - ANGLE;
    sign = 1;
}
else if((ANGLE >= SINE_TABLE_LEN2) && (ANGLE < SINE_TABLE_LEN3))
{
     index = ANGLE - SINE_TABLE_LEN2;
     sign = -1;
}
else if((ANGLE >= SINE_TABLE_LEN3) && (ANGLE < SINE_TABLE_LEN4))
{
     index = SINE_TABLE_LEN4 - ANGLE;
     sign = -1;
}
else
{
     ANGLE -= SINE_TABLE_LEN4;
     index = ANGLE;
     sign = 1;
}
audio_sample_l = sign * sine_lut[index-1];
audio_sample_r = audio_sample_l;

I initialize my variable ANGLE at 0 and FREQUENCY_VALUE at 1. When a plot the fft of my sinus I get some odd harmonic noise that I can not explain. When I compare the output samples with the one from my lut, they are the same (no index error). The problem may come from the fact that I do not start at the beginning of my lut but at the previous sample to avoid phase changes. There may be some rounded error also. My third lead is the initialization values and the index calculation (test in the loop also). What do you think ?

Anyone have some ideas to change that ?

Trace your logic when ANGLE is incremented to SINE_TABLE_LEN4 -- you'll find you are accessing below the buffer.  You should really adopt a zero-based index scheme when coding in C.


how can we align wto signals with a time lag in between? I have two signals from two different sensors. I cropped the area i am interested in from both sensor signals. The problem is that the signals start at different time instants. I am searching for a way to align the signals just like below. It is technically just moving the blue signal to align with orange one. I tried using cross correlation but it gave a weird result. Any help please.
 first image:  this is the result of cross correlation of two signals
Second image :after applying code on cross correlation to align

I thought you might have to smooth, sort or create an envelope of your signal, which might help. But then i found this. This should help you.

Or if you wish to use MATLAB, you can use this command.


Practical examples of ARMA model I am studying the Kalman filter and its basic implementation, and it was asked to use the filter to estimate a signal observed in noise
$$y(n) = x(n) + v(n)$$
where $v(n) \sim \mathcal{N}(0, \sigma^2)$ and $x(n)$ is modeled as an ARMA(2,2) process
$$x(n) = b_0 u(n) + b_1 u(n-1) - a_1 x(n-1) - a_2 x(n-2).$$

The exercise is conceptually intersting, but I was wondering if there is an direct application for this case. I mean, which processes, in real life, can be modeled as an ARMA(2,2) and in which context the Kalman filter is used, not as a predictor, for a signal like $s(n)$?

Thanks!

ARMA models are useful when you need to model a Signal plus Noise situation where the signal is an AR process and the noise models sensor noise. The overall model is an ARMA model.

See HL Van Trees, Detection, Estimation, and Modulation Theory, vol 4 Array Processing. He gives an example of a Spatial AR process sensed by noisy sensors. The overall model is ARMA

ARMA model are also useful in situations where you have strong nulls n your spectrum like multi path and you have a noise like signal


What are poles and zeroes (with respect to the inputs and outputs of a system)? I get it, about poles and zeroes, when we talk in terms of the transfer function.

But, if the transfer function is the ratio of the output to the input, then if the input signal is zero, then the gain should go to infinity, right?

Am I missing something?

A transfer function is defined as the ratio of the transform of output to the transform of input where all initial conditions are zero.


DFT frequency resolution exercise I have a discrete signal that goes as follows:
$$x[n]=[-1,4,-1,0]$$
I have already done the DFT for the signal, with the following result:
$$X[0] = 2, X[1]=-4i,X[2]=-6,X[3]=4i$$
But for some reason, I seem to have trouble getting the frequency resolution for these samples(?).
My sampling interval (T) is 0.03125 seconds, and with some logical deducting I came up with these following frequencies:
$$f_0=0,f_1=8,f_2=16,f_3=8$$
However, I do not understand how I came to this conclusion, but it seems correct. I tried to use the following formula:
$$\Delta f = \frac{1}{T_{0}} = \frac{1}{NT}$$
Could someone explain how I ended up with these frequencies?

Carry the units and it becomes understandable.

T =  0.03125 seconds per sample

Thus your sampling rate is

1/T = 32 samples per second

Each bin corresponds to the cycles per frame.  Your frame has four samples.  Suppose the bin index is k.

k (cycles per frame) * 32 (samples per second) / 4 (samples per frame) = 8k cycles per second = 8k Hz

So the frequencies for your bins, in Hz, become

k   f_k
0   0
1   8
2   16
3   24

However, the DFT wraps around, and choice of the principal frequency among many aliases is a matter of convention.  It is mathematically sounder to consider the upper half of the DFT as negative frequencies.  For real valued signals, the upper half is just the complex conjugate of the lower half and thus provides no new information so it is usually ignored.  So, if you want to consider f_3 as 24 or -8, it doesn't really matter.  For a real valued signal, a negative frequency is the same as a positive frequency.


How Is the MATLAB `filter()` Function So Fast? How is the filter() implemented in MATLAB? How is it so fast?
What's the fastest implementation of a FIR filter?

Actually it seems MATLAB implementation of the filter() function is pretty straight forward and not fast.

For a fast implementation, have a look at FilterM by Jan Simon.

Update

In the latest releases of MATLAB (From R2016b and above) the performance of the filter() function has improved.
The metdhos to accelerate those operations are usually based on:

Utilizing SIMD Operations (SSE / AVX in x86).
Multi Threading.
Loop Unrolling.


Homography for (small baseline) Stereo Image Registration If the distance between two stereo cameras (baseline) is much smaller than the distance to the geometry, is it possible to align the resulting two stereo images using a homography?

For most cases (and I tend to say) yes, but this also depends on the geometry. If there is a perpendicular plane and two cameras observe different sides, the assumption will cease to hold.
Also, in the cases where objects are far and no crazy things are going on, you could also parameterize the intrinsics using an affine camera model.


Example exercise - Shannon's law on PCM-digitized signal I am studying for a Digital Communications course in my university, and I am presented with the following example exercise.

  You have a sound signal with spectral components in the 300-3000 Hz range. The signal is going to be digitised using Pulse Code Modulation.

  a. For SNRdb == 30dB, what is the number of uniform quantum levels needed?

  b. What is the minimum data transfer rate needed for the transmission of the digitized-with-PCM signal?

I know that for both of those questions I need to use Shannon's law, C = B*log_2(1+SNR) with SNRdb = 10*log10(SNR).

For task a. I need to find M, for task b. I need to find C.

My question: What do I use as B in both of the questions? How do I utilize the digitized-with-PCM info?

Do I take a median 1650Hz in a?

In b, do I write that PCM uses from 8 to 20something bits per sample, and use B = 20*3000Hz? Or do I just use the maximum B of the range, 3000?

Any help would be very much appreciated.

I know that for both of those questions I need to use Shannon's law, C = B*log_2(1+SNR) with SNRdb = 10*log10(SNR).

I'm not quite sure that's true. For the first question,

  For SNRdb == 30dB, what is the number of uniform quantum levels needed?

quantization theory is pretty sufficient: In every digitally represented signal, there's so called quantization noise; it's nothing but the power in the difference between analog and quantized signal. For example, every additional bit in the ADC halves the size of the quantization step, reducing the average voltage error by one half, leading to four times less noise power (or, equivalently, 6dB better quantization noise).

So, this question, if you ask me, is: "When you need to achieve quantization noise better than 30dB, how large can your quantization levels be?".

Solution is relatively simple: Let $P_s$ be your signal power. Your Quantization-limited SNR is then

\begin{array}
\mathrm{SNR}_Q &= \frac{P_s}{P_\text{quantization noise}}
\end{array}

Without loss of generality, we set $P_s=1$. Hence, with $q$ the number of quantization steps and 1 the overall ADC range,

\begin{array}\\
{\mathrm{SNR}_Q}^{-1} &= {P_\text{quantization noise}}\\
&= E\left[\left(x_\text{digital} - x_\text{real}\right)^2\right]\\
&= \text{var}({e_Q})\text{, }e_Q\text{ being the quantization noise Random Variable (RV),}\\
&= \int\limits_{-\frac1{2q}}^\frac1{2q}{\left({e_Q}^2q \right)d{e_Q}}\text{, since variance of a uniform RV over }[-\frac1{2q}; \frac1{2q}]\text{ is known,}\\
&= q \left[ \frac13 {e_Q}^3 \right]_{-\frac1{2q}}^\frac1{2q}+C\\
&= \frac q{3}\left(\left(\frac1{2q}\right)^3 -\left(-\frac1{2q}\right)^3  \right)+0\text{, }C\text{ being zero because of plausability}\\
&= \frac q3\left(\frac1{8q^3}+\frac1{8q^3}\right)\\
&= \frac q3\frac1{4q^3}\\
&= \frac1{12} q^{-2}\\
&\overset!\le {30\text{dB}}^{-1}\\
&=\frac1{1000}\\
&\iff\\
q^2&\ge \frac{1000}{12}\\
&\iff\\
q&\ge \sqrt{\frac{500}6}\\
q&\ge 9.1\text{something}\\
q&\ge \left\lceil{9.1}\right\rceil \\
&=10
\end{array}

assuming that the ADC is perfectly linear and thus has absolutely uniform quantization noise for a signal that has uniformly distributed power over the whole ADC range, and thus, the error is uniformly distributed over the quantization step of size $\frac1q$.

That was fun.

  b. What is the minimum data transfer rate needed for the transmission of the digitized-with-PCM signal?

Ha, Shannon, here we come.

PCM means we can't trick around with undersampling. So, our signal's bandwidth is 3kHz, meaning we need to sample with 6kHz. Every sample has information $I=-\log_2(P(x))$, $P(x)$ being the probability that value x is represented in the sample. Still assuming the signal is uniformly distributed (the assignment did not say it was a sum of sinusoids or similar!), $P$ is constant, $P=\frac1q$, hence $H(X)[\text{bit}]= -\log_2{\frac1{10}}= \frac{-\log_{10}{\frac1{10}}}{\log_{10} 2}=\frac1{0.3\dots}\approx 3.32\text b$.

So the bandwidth (bitrate) we need to transport is

\begin{align}
r_b &= 3.32\text b \cdot 6 \text{kHz}\\
&= 19.93\frac{\text{kb}}{\text s}\text{ .}
\end{align}

So with approximately 20kb/s, you can do it.

Those are very good questions. I'd like to ask you to post your own answer if you get an official solution, of course only after getting an O.K. from the professor/lecturer/TA to share his intellectual property (maybe pointing him to this would help). (I of course still like if you accept and upvote this one. It could also certainly use some improvement/corrections!)


Implementation of gain scheduling for a PI controller I'm designing a PI controller for a boost converter. I was able to find a set of gains that fit for all situations, but I get a slow response at full-load.

I was thinking instead of using gain scheduling, 1 set of gains for little to no load and 1 set of gains when load is significant (say more than 25% of the maximum load).

Now my question. What are the gotchas about gain scheduling? In books they don't say much about implementation details.

I was thinking that in order for gain scheduling to work, I need to apply the Ki gain of the integrator AFTER the integration and not before.  That way I canchange seamlessly from one set of gain to the other.

I.e :

Use this =>
integrator = integrator + error;
output = Kiintegrator + Kperror;

Instead of
integrator = integrator + Kierror;
output = integrator + Kperror;

Am I on the right path?

Regards

Benjamin

if you want the transition between two sets of controller parameters to be smooth, you need to apply the integral gain before the integrator. The controller output will perform a step otherwise, which is the exact opposite of seamless.

Your speaking of a PI controller and you have to apply a suitable strategy for the proportional gain too. Possible solutions are, but not limited to:

apply a 2nd order filter to the controller parameters to ensure smooth transistion, this might affect controller stability during transistion but this shouldn't be a problem with your simple control loop.
Implement a controller per parameter set, one of the controllers is active while the others are in tracking mode, switch over when another set of parameters is required. This will ensure, that the controller runs instantly with the correct parameters which typically doesn't destabilise the control loop. (Simulink PID Block, search for tracking)
Try to find a function that expresses the perfect controller parameters based on the current load situation and adjust the parameters accordingly. This might be tedious though.


What interpolation methods can I use to get the tightest fit for these curves? I am working with MRI images of the brain that have certain areas marked by hand like  and .

I am trying to come up with an interpolating function that will let me describe such curves so that I can use machine learning techniques to automatically generate such markings on unmarked images like  and .

I've considered using cubic spline interpolations to approximate the curve, but as far as I know, I would have to split the curve into multiple parts to interpolate it.

Am I on the right track or are there other (better?) ways to do this?

If you have a time series of images, I recommend the CONDENSATION algorithm for the interpolation of curves (see also the wikipedia page). There are a number of implementation available, notably in opencv.


How does ICA handle inevitable delays in signals? I am currently reading up on and teaching myself ICA from a number of good sources. (Also see this post for past context). I have the basic jist down, but there is something I am not clear about.

For a scenario where multiple signals are impinging on multiple spatial sensors, (of course, with number of sensors >= number of signals), it is inevitable that for any one sensor, all the signals arriving to it will have different delays/phase-offsets associated with them, as compared to those arriving at a different sensor.

Now, as far as I know, the signal model for ICA is a simple mixing matrix, where the total energy arriving at any one sensor is modeled as nothing but a simple linear combination of all the other signals of interest. Every sensor has a different array of linear combination coefficients associated with it. So far so good.

What I do not understand, is that inevitably there are going to in fact be some delay/phase-offset among individual signals arriving at individual sensors that differ from one another. That is, $s_1(n)$ might arrive at $sensor_1$ at some time 0s, while that same $s_1(n)$ arrives at $sensor_2$ attenuated, but also at some delay or phase difference. The way I see it this is physically inevitable.

...How can it be that this is not modeled in the mixing matrix? It seems that delays will make a huge difference. No longer are we talking about simple linear combinations anymore. How does ICA handle this? Have I missed something here?

I should also add as an addendum, if indeed ICA cannot handle delays, then what applications does it find usefulness in? Clearly spatial ones with sensors are out!

Thanks

One of the most successful uses of ICA has been in the study of electrophysiology (i.e. brain activity), principally EEG (Electroencephalography) and MEG (Magnetoencephalography). They are used to remove artefacts (such as electrical impulses caused by muscle movements (eye blinks etc)) without the need for reference channels. In this application the spatial separations between sensors is minute compared with the propagation speed of the waves, and as such the assumptions of ICA effectively hold.

For fMRI, which relies on blood flow in the brain, the temporal delay issue is more significant. One approach, taken in the paper Latency (in)sensitive ICA. Group independent component analysis of fMRI data in the temporal frequency domain by Calhoun et al (2003) attempted to solve this problem by making estimates of the time delay in each voxel, and then using this as prior information in a modified ICA. Maybe something like this could be applied in your domain?


FPS to Hertz in a camera Not sure if this is the proper place to ask this, but I'm sure someone here will know the answer.

I have this optical camera (TrackIR 4), which, for what you can see on the table, samples images at a rate of 120 frames per second.

Now how would you do a conversion from fps to Hertz to find out the frequency rate at which the camera samples data?

Maybe it's a stupid question, but I wanna make sure I get the right value.

The physical definition of frequency is the number of oscillations per second and measured in Hertz (Hz) or more general: 1/(duration of one period of a periodic event in seconds). In your case the periodic event is a frame and it occurs 120 times per second so the frequency is 120 Hz.

If you want to find out the sampling frequency of the ADC device in the camera you could calculate it like that:
$$
f_\mathrm{s} = N_\mathrm{pixel} \cdot R_\mathrm{frame}
$$
where $N_\mathrm{pixel}$ is the number of pixels in the sensor and $R_\mathrm{frame}$ is the frame rate. With the numbers from the data sheet I obtain a sampling frequency of 12.3 MHz. I can't tell if this is the actual frequency used by the camera. If it uses oversampling or interlacing the value could differ.


Spatial and temporal frequency of a wave I've sampled a propagating wave so that I know its amplitude at several equidistant points on a line, and at several points in time. That is, I have this discrete function:

f[t,x]

Now, I want to find the frequency of the wave. But there seems to be two frequencies I can calculate: by finding the FFT at different points in time (e.g. FFT{f[1,x]} ) and at different points in space (e.g. FFT{f[t,1]}).

Only the latter ("temporal frequency") gives me a main frequency in Hz. The other one gives one in 1/m, but I suppose that if I knew the wave propagation speed in m/s, I could convert one to the other. Is this true?

All of the above I've just reasoned through myself. Could someone point me to a source where I could read about the two different types of frequency?

Edit: I'll add that the time-space relation is "normal", i.e. the wave propagates like $g(kx-\omega t)$. Hope I'm not being too unclear. I realise that I'm lacking the terminology needed to explain what I mean.

Edit 2: It might be important to know that $f$ is a pulse that is localised in time and space.

Edit 3: I'll put this question here too: Is there a way to take advantage of all the sample points (a 2D matrix) to estimate the central frequency, assuming I know the wave speed?

The two frequencies you are referring to are the spatial frequency and temporal frequency of the wave, and you are correct in your reasoning on converting one to the other. The spatial frequency refers to how many complete periods the signal goes through for a given unit of distance (eg. cylcles/m) while the temporal frequency refers to how many complete periods the signal goes through for a given unit of time (eg. cyles/s or Hz).

The physical meaning of the two frequencies should be fairly easy to understand if you imagine the wave you are sampling is a simple sine wave. Then it is easy to see that for any constant point in time, your signal will be a frozen sin wave which is periodic in space, and thus has a "spatial frequency". Conversely, for any constant point in space, the value of the signal will oscillate over a regular period in time, and thus have a "temporal frequency" component. You can also imagine that the faster the wave was travelling through space, the faster the signal would change at any one point in space, giving rise to a higher temporal frequency and underlining the relationship between the two properties via the wave's speed.

Regarding using all the sample points to estimate the wave frequency, I believe you should be able to use the known speed to scale one of your dimensions to match the other, since $distance=velocity*time$. Then a 2D FFT of the matrix should give you information based on all of the available samples in whichever dimension you choose.

This article gives a brief description of the relationship between the two in the context of radar imaging and this article talks about the use of the 2D FFT for analyzing higher dimensional signals.


How to graph magnitude plots for basic 4-pole filters? I have developed magnitude plots for basic one-pole and two-pole filters derived from various Physics PDF's and tutorials I have found online and calculated from:

https://www.desmos.com/calculator/mjudrnexbo

But I cannot figure out what the correct plots for any simple four-pole resonant LPF, HPF, BP, & BR filters would be in the same terms.

Any help on what equivalent simple four equations would be? Thanks.

For Butterworth LP and HP

https://www.desmos.com/calculator/dxcnvc63av

Play the variable n.

Sources:
https://www.electronics-tutorials.ws/filter/filter_8.html
http://fourier.eng.hmc.edu/e84/lectures/ActiveFilters/node6.html


Alternative to Orthogonal Matching Pursuit (OMP) Algorithm In the Compressed Sensing context, assume there is a signal $ x \in {\mathbb{R}}^{n} $ which is $ k $ sparse. Namely its Pseudo $ {\ell}_{0} $ Norm is $ {\left\| x \right\|}_{0} = k $ (The signal has only $k $ non vanishing elements) where $  k << n $.

Given a Model Matrix $ A \in {\mathbb{R}}^{m \times n} $ the measurements are given by (This is the Model):

$$ y = A x $$

The recovery problem is given by:

$$ \arg \min_{x} {\left\| A x - y \right\|}_{2}^{2} \; \text{s. t.} \; {\left\| x \right\|}_{0} = k $$

Since the exact solution to the problem above is hard to find the recovery (Estimation) of the signal $ x $ from the measurements $ y $ is usually done using Orthogonal Matching Pursuit (OMP) Algorithm.

Basically the OMP finds iteratively the elements with highest correlation to the model.

The question is, since the measure of the quality to select indices is based on Correlation, why can't one just find the support using the following:

[vCorrVal vCorrIdx] = sort(A' * Y);
vSignalSupport = vCorrIdx(1:k);

The result is almost the same as that of normal OMP.

The main advantage of OMP is that the residual is orthogonal to the current solution.

Let's say you select all $k$ columns from $A$ (also called atoms) at once and let us also presume that $A$ is an overcomplete basis (this is more or less the standard in OMP literature).

Now, with your method, if the atom that correlates the most with your measurements $y$ is linear-dependent with $p < k$ other atoms in $A$ you will end-up with an $k-p$-sparse signal, because $p$ entries will be more or less redundant. The same argument can of course be extended to less correlated atoms. You might also be lucky and never see the phenomenon.

Let's take the same example but with OMP this time. During the first iteration you would select the atom that correlates the most with measurements $y$. After that you compute the coefficient in $x$ such that the new residual is orthogonal to the current measurements approximation. In other words you got the most of the information provided by the currently selected atom so during the next iteration you are very likely to pick an atom that contains fresh information (ask yourself what would happen with linear dependent atoms in this case).

Here is a list of atom selection look-ahead strategies based on OMP and OLS that you might find interesting to read: POMP, LAOLS and POLS.


Decimate signal by 64 I have a 16bit SAR ADC sampling at 6.4kHz. The Anti-Aliasing Filter is a active 4th order Bessel lowpass (2x LTC6081) with -3dB at 200Hz and -80dB at 3.2kHz.

Now I want to decimate these 6.4kHz to 100Hz. The uC is a 8bit RISC (Atmel AVR) running at 20Mhz.

My first idea was multiple digital bessel biquad filters and downsampling by a factor of 2 but this performs horribly. Then I read about half-band and polyphase filters but I don't have any experience in these topics.

So my primary question is: how should I downsample my 6.4kHz signal to 100Hz?

So, yes, since your decimation is a power of two, going through cascaded half-bands is a clever thing to do!

First of all, what's special about a half-band filter? If you implement a filter of 1/2 bandwidth as FIR, then you can choose to design it so that the transition width's center is exactly at ¼ the sampling rate – that allows for symmetry, and that in turn gives you a FIR filter where every other coefficient is zero.

Zero coefficients don't need to be summed. This gives you the steepness of a $K$-tap filter at the cost of $\frac K2$!

Now, let's look at the fact that you're decimating by a factor of 2:

That means that you throw away every other output sample. Polyphase decomposition of your filter allows you to write your single filter as a sum of two alternatively "fed" filters, so that instead of calculating $K$ products and sums per time step, you calculate $\frac K2$. That's another effort reduction by a factor of 2.

So, what you do is:

you pick an existing decimating FIR filter implementation in C or assembly, or write it yourself (it's but a simple for-loop). GNU Radio has one; it's in C++, and uses optimized filter kernels, but it might help illustrate how you split up one filter into $K$ subfilters to reduce the load by a factor of $K$ in a $K$-decimating FIR.
You design a half-band filter. You probably don't need something fantastic; make sure you get as many zeros in there as sensible.
you concatenate $\log_2 64=6$ of these (in fact, that's only nearly efficient; read that article, it's good).

General comment: Using an 8-bit microcontroller to process 16 bit numbers is a bad idea through and through. Atmegas are usually more expensive than similarly peripheral-equipped ARM counterparts, and really, they aren't even remotely adapted to DSP loads. Don't do that to yourself. This task is easy even on a cheap 20 MHz ARM Cortex-M0, and you could have a Cortex-M4F for the price of a larger ATMega. That would have a floating point unit; not that you need that for any of this (you can do all these calculations in fixed point, like you'd have to do it on your AVR), but it does make one's life easier further down the processing chain. I assume you use an AVR "because you already have that, and have experience with that"; but believe me, you're shooting yourself in the knee. More suitable MCUs are cheap, and so are development boards (they literally start at 2$) sufficient to connect your ADC to them.


How do I calculate the step response of a discrete-time system? I am given the impulse response

$$h[n]=2\left (\frac{1}{2}  \right )^n\sigma [n]-\delta [n]$$

where $\sigma$ is the unit step function and $\delta$ is the unit impulse function. I know that

$$a[n]=\sum_{k=-\infty}^{\infty}h[k]\sigma [n-k]=\sum_{k=-\infty}^{n}h[k]$$

Thus, my step response should be

$$a[n]=\sum_{k=-\infty}^{n} \left( 2 \left (\frac{1}{2}  \right )^k \sigma [k]-\delta [k] \right)$$

How do I calculate this?

Taking the $\mathcal Z$ transform of the impulse response,

$$H (z) = \dfrac{2}{1 - \frac 12 z^{-1}} - 1$$

Hence, the $\mathcal Z$ transform of the step response is

$$S (z) := \dfrac{H (z)}{1 - z^{-1}} = \cdots = \dfrac{3}{1 - z^{-1}} - \dfrac{2}{1 - \frac 12 z^{-1}}$$

Taking the inverse $\mathcal Z$ transform, we obtain the step response

$$s (n) = \begin{cases} 3 - 2 \left( \frac 12 \right)^n & \text{if } n \geq 0\\ 0 & \text{if } n < 0\end{cases}$$


How is a realtime BiQuad/IIR possible? I was thinking, if a realtime system is one in which the processing time is shorter than the sampling period, how does a feedback filter such as a biquad perform the computation without waiting for, say, in the case of the biquad, two additional samples?

You are confusing "processing time" with "latency".

Real-time filters are able to generate output samples at the same rate as they receive inputs.
They would however induce latency, meaning that the bulk of the energy generated for an input appears later in time with respect to that input.

Consider for example the following input/output:

time  : 0 1 2 3 4 5 6 7
-----------------------
input : 1 2 3 2 1 0 0 0
output: 0 0 1 2 3 2 1 0

The output is generated at the same rate as the input is received, so the system can perform in real-time. However, there is a latency of 2 samples in the output when comparing with the input.

Whenever the expression of the filter relies on past samples which have not yet been received (at the start of an input sequence), initial conditions are assumed (often assuming that those previous values were zeros).


The need for phase wrapping? I understand the need for phase unwrapping and the complexity behind it. But, I'm struggling to understand why the phase information gets wrapped between $-\pi$ to $\pi$ in the first place?

I'm not sure if I understand your problem, but I'll give it a try. If you have a complex frequency response

$$H(\omega)=|H(\omega)|e^{j\phi(\omega)}\tag{1}$$

where $\phi(\omega)$ is the phase, then you probably know that for a given frequency $\omega_0$, the frequency response might as well be written as

$$H(\omega_0)=|H(\omega_0)|e^{j(\phi(\omega_0)+2k\pi)}\tag{2}$$

where $k$ is some integer. The question is, when you compute the complex value $H(\omega_0)$, which of the infinitely many possible phases do you choose? The convention is to choose the principal value which is in the interval $(-\pi,\pi]$. The consequence is that whenever you compute the function $\phi(\omega)$ for continuous $\omega$, you'll get jumps as soon as the phase hits $-\pi$ or $\pi$, but no matter how you define the phase, you will not be able to avoid jumps. So if you want a continuous function you need to do phase unwrapping.


Understanding this Haar Wavelet Example I'm trying to understand Wavelet Transform , I've done well enough to understand Continuous Wavelet Transform which was easy enough, where we simply stretch the wavelet and match it with the original signal.
But i couldn't be able to understand the Discrete Wavelet transform in fact I've always had a trouble understanding Discrete of any transform (like DFT , DCT which are still not clear to me).

According to Wavelet transform we Stretch our signal and then match it with the original signal to find out the frequency.
This is the example shown in the book,

  In this hypothetical example the student does fairly well the first
  half of the term then neglects his or her studies for the last half.
  Thus the exam scores for the term were 80%, 80%, 80%, 80%, 0%, 0%, 0%,
  and 0%* We can tell the average of all the scores (40%) and when the
  scores “tanked” after the 4th exam just by looking. Knowing the answer
  in advance, however, is a good way to learn and to verify the wavelet
  transforms. Then we can use them with confidence on real-world data
  where we can’t simply “eyeball” the final values.

  We will now walk through the CWT process step by step using the
  simplest of the wavelet filters on this example. We begin by comparing
  the humble Haar wavelet filter, [1 –1]

Signal -> [80 80 80 80 0 0 0 0]
Filter(or signal that will be stretched later) -> [1 -1]

Comparing the first 2 points with the wavelet filter we obtain 80 – 80 = 0. For
this very simple high-pass filter we can say there was no change in the first 2
exam scores.

This is how we've done it,

80*1 + 80*(-1) => 80 - 80 = 0

Now We have,

[0 80 80 80 0 0 0 0]

Shifting once to the right and applying filter again we get,

80*1 + 80*(-1) => 80 - 80 = 0

Now We have,

[0 0 80 80 0 0 0 0]

Shifting once again the right and applying filter again we get,

80*1 + 80*(-1) => 80 - 80 = 0

Now We have,

[0 0 0 80 0 0 0 0]

Shifting again,

80*1 + 0*(-1) => 80 - 0 = 80

Now We have,

[0 0 0 80 0 0 0 0]

We will eventually end-up with

[0, 0, 0, 80, 0, 0, 0]

This is significant in that this wavelet process of comparison and
shifting has just indicated a large change between the 4th and 5th exam. We
have “found the discontinuity”.

This actually makes sense as we had a large change from 4th to 5th point But how does this discontinuity defines a frequency here ?

Now if we stretch our signal from [1 -1] to [1 0 -1] and repeat the process,

We will eventually end-up with,

= –80 –80 0 0 80 80 0 0 0 0

When we used filter [1 -1] , our output at least makes some sense it clearly shows a rapid change, but what happened here in case of 3 point filter ? what are these weird values (-80 -80) ?

If we keep scaling to 10th times here is what we'll eventually get,

(scale = 10) –320 –160 0 160 320 320 240 160
(scale = 9)  –240 –80 80 240 320 240 160 80
(scale = 8)  –320 –160 0 160 320 240 160 80
(scale = 7)  –240 –80 80 240 240 160 80 0
(scale = 6)  –240 –160 0 160 240 160 80 0
(scale = 5)  –160 –80 80 160 160 80 0 0
(scale = 4)  –160 –80 0 80 160 80 0 0
(scale = 3)  –80 0 0 80 80 0 0 0
(scale = 2)  –80 0 0 0 80 0 0 0
(scale = 1)   0 0 0 0 0 0 0 0

While studying the Continuous Wavelet Transform the scaling made sense but i could't understand the scaling here ? what are these points -80, -160, -240, -320 represents ?

This is the final output,

How can i relate this output from my original signal ?

It looks like there is a problem with your scaling.

The scaling for the DWT has the same interpretation as it does for the CWT.

For the DWT, the scale of the analyzing function (wavelet) is increased using a dyadic scale (increasing by factors of 2) in non-overlapping time intervals.  Increasing scale corresponds to narrower frequency distributions in the frequency domain.  Basically think of the analyzing function (wavelet) as a bandpass filter for which you are both narrowing the bandwidth and shifting down the center frequency by increasing the time scale.

As an example (see images below) I've transformed an audio signal using a CWT and then again using a DWT.  Notice that the results look similar, they scale in the same way, except in the DWT, the regions of analysis are discrete non-overlapping intervals (sharp boundaries).  Notice also that as the time scale increases (heading south on the CWT diagram), you have fewer DWT coefficients, so the number of blocks in the image decrease.  I don't see this kind of scaling in your data.


Calculating the true values of biquad filter coefficients from normalized form I have a C code in which the IIR filter(BiQuad) uses only 4 (normalized) coefficients per stage (a2/a0,a1/a0,b2/b0,b1/b0) instead of 6 coefficients (a0,a1,a2,b0,b1,b2). I know that the value of a0 is always 1. So, it is easy to find the true values of the a2,a1 coefficients (since a0 value is 1).  But, i wanted to know the true values of b0, b1, b2.

So, I just want to know is there any way to retrieve/derive the true values of "b" coefficients?  One of my friends told me that we can get the values by using Matlab / octave? But, I have no idea how to do it?

Note that a biquad has $5$ degrees of freedom (not $6$), because $a_0$ can always be chosen as $a_0=1$ without loss of generality:

$$\begin{align}H(z)&=\frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}\\&=b_0\cdot\frac{1+\hat{b}_1z^{-1}+\hat{b}_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}},\qquad \hat{b}_1=\frac{b_1}{b_0},\;\;\hat{b}_2=\frac{b_2}{b_0}\tag{1}\end{align}$$

Eq. $(1)$ shows that $b_0$ is just an overall gain (for fixed values of $\hat{b}_1$ and $\hat{b}_2$), as you've noted yourself.

If you don't know the formula that produced the filter coefficients, there is generally no way to know the intended gain. As mentioned in Dan's answer, the gain is often chosen such that quantization effects are minimized. It can often be compensated for after all computations have been performed.

In some cases, however, it is possible to estimate the intended gain. E.g., if the biquad is a low pass filter, it is reasonable to assume that the original frequency response has a value of $1$ at DC. This means that

$$b_0\frac{1+\frac{b_1}{b_0}+\frac{b_2}{b_0}}{1+a_1+a_2}=1\tag{2}$$

Similarly, for a high pass filter a common scaling makes sure that the frequency response equals $1$ at Nyquist, which implies

$$b_0\frac{1-\frac{b_1}{b_0}+\frac{b_2}{b_0}}{1-a_1+a_2}=1\tag{3}$$

You can come up with similar guesses for the intended scaling for other standard filter types, like band pass filters (unity gain at the center frequency) and band stop filters (unity gain at DC and/or Nyquist).


References ISAR - Inverse Synthetic Aperture Radar Anybody know a reference to start understanding ISAR? Book or article.

Thanks!

A couple of books I'd recommend that are similar to the Richard's book

Stimson's "Introduction to Airborne Radar" has a chapter on ISAR.
For a bit more detail you can look at
Wehner's "High Resolution Radar"

If you look at the Polar Format algorithm in SAR (I'd recommed Soumekh's - Synthetic Aperture Radar Signal Processing) it is very similar to the ISAR algorithms after you corrected the target's translational motion (this is not necessarily trivial). Soumekh's book doesn't really address the ideas of motion compensation, target tracking, and range cell migration - these are essential concepts if you are implementing an ISAR algorithm on Airborne data. If you are just trying to get an idea of ISAR processing then you just need a little bit about them.

A good book that covers a wide range of ISAR topics is:"Inverse Synthetic Aperture Radar Imaging: Principles, Algorithms and Applications" by Chen and Martorella.

One piece of advice, if you are reading papers - pay very close attention to the type of transmitted signal and the type of range compression that is used. Often authors just state the signal after "preliminary" processing is given by ... . If you do deramp on receive vs traditional matched filter, the nature of the signal is very different. Similarly, if a stepped frequency waveform is used then the resulting signal processing is different. This is one of the reason's I liked Soumekh's book.

One other book is:Inverse Synthetic Aperture Radar Imaging with MATLAB Algorithms by Caner Ozdemir. This book explains the concepts well, but I was unable to follow the mathematical development and the Matlab code is terrible. You can't easily modify it change the processing parameters.


Delayed impulse does not delay spectral response Inspired by Dr. Sarwate's answer to the question  How replicas are formed in Frequency domain when a signal is sampled in Time Domain? , I decided to see for myself the underlying concept.

I interpreted Dr. Sarwate's answer that sampling was a sequential activity. All one needed to do was to multiply an impulse with the signal meant to be sampled. And this needed to be repeated with delayed/shifted versions of the impulse.

Multiplication in time domain is convolution in frequency domain and vice-versa. I thought that with the impulse shifting by every sample, its spectral response would also undergo a similar shift in frequency domain. But it seems not to be the case.

I wrote a toy program in octave and the plots of the frequency response havent brought any joy. What have I misunderstood here?

## Num Samples
N = 10;

## Row vector representing samples (all zeroed out)
x = zeros(1,N);
y = zeros(1,N);
z = zeros(1,N);

## Undelayed impulse (For lack of a better word)
x(1) = 10;

## Delayed by 1 sample
y(2) = 10;

## Delayed by 2 samples
z(3) = 10;

## Obtain spectrum of undelayed impulse
s_spectrum_x = fft(x);
s_spectrum_abs_x = abs(s_spectrum_x);

## Obtain spectrum of impulse delayed by 1 sample
s_spectrum_y = fft(y);
s_spectrum_abs_y = abs(s_spectrum_y);

## Obtain spectrum of impulse delayed by 2 samples
s_spectrum_z = fft(z);
s_spectrum_abs_z = abs(s_spectrum_z);

>>disp(s_spectrum_abs_x);
10 10 10 10 10 10 10 10 10 10

>>disp(s_spectrum_abs_y);
10 10 10 10 10 10 10 10 10 10

>>disp(s_spectrum_abs_z);
10 10 10 10 10 10 10 10 10 10

Shouldn't the initial bins of y and z be 0?

The convolution theorem ("Multiplication in time is Convolution in frequency") states:

$$\mathcal{F}\{x\cdot y\}=\mathcal{F}\{x\}*\mathcal{F}\{y\}$$

where $x, y$ are two time-domain signals, $\cdot$ denotes element-wise multiplication and $*$ is the convolution operator. So, the convolution in the frequency domain, is carried out between the Fourier Transforms of the original signals.

As an example, let us use your impulses with frequency response:

$$\mathcal{F}\{\delta(t)\}=1$$

Further, note the time-shifting property of the FT:

$$\mathcal{F}\{x(t-t_0)\}=\mathcal{F}\{x(t)\}\exp(-j2\pi ft_0)$$

where $x(t)$ is any signal you like.
Now, you can readily explain the results you obtain for your shifted impulses:

$$\mathcal{F}\{\delta(t-1)\}=\mathcal{F}\{\delta(t))\}\exp(j2\pi f)$$

where $f=[0/N, 1/N, 2/N, \ldots, (N-1)/N]$ (due to the discrete FT).  Similar for the two-shifted case:

$$\mathcal{F}\{\delta(t-2)\}=\mathcal{F}\{\delta(t)\}\exp(-j2\pi f2)$$

Note that $\mathcal{F}\{\delta(t)\}=1$. You display the absolute value of your spectra: $\left|X(f)\exp(j2\pi ft_0)\right|=|X(f)|$ for any $X(f)$ you like. That's why you see three times the same output. Try plotting real and imaginary part separately, then you will see the effect of time-shifting. You might also have a look here which is an article I wrote about the effect of the time-shift and frequency-shift.


detection of periodicities in n-dimensional signals Generally speaking, what analyses are necessary and sufficient for the detection of periodicities in an n-dimensional signal amounting to a discretely sampled density distribution over n-dimensional space?

Assume that the signal-to-noise ratio and sampling frequency in all dimensions are adequate for detection of any periodicities that are present.

Would an n-dimensional Fourier transform and an n-dimensional autocorrelogram suffice?
Would these be limited to detecting periodicities that are aligned
with one of the basis vectors of the n-D space?

Suggested references
   for further reading would be appreciated.

Would these be limited to detecting periodicities that are aligned with one of the basis vectors of the n-D space?

Well, if you say "sampling is sufficient", then it follows that the whole space is the span of the base vectors (in fact, you usually have them orthogonal).

From that follows that a non-base-vector aligned oscillation can be represented by a linear combination of oscillations that are aligned with base vectors.

From that follows that yes, with your requirement that sampling frequency is sufficient to begin with, the usual methods of 1D frequency analysis work.

Note that you can't generally detect all periodicities with just a sufficiently fast sampling – if something is periodic with a period  longer than you observation window, you have no information about that. So, observation length must be sufficient, too, unless your signal model restricts both highest and lowest frequency.


How do spectrum analyzers detect peaks? Using peak search, I'm able to put the cursor on any of the several peaks on the spectrum analyzer display.  I am trying to do something similar in software, with the output of the FFT of the radio spectrum.  Using normal peak detect functions (such as those included in Scipy) does not seem to work.  I think that there must be a particular algorithm used for radio, where particular circumstances exist (relatively flat noise floor, etc.) , what is it?

EDIT:

Here is a screenshot:

The orange markers are where the peaks are thought to be.

How I made this was filtering and using scipy.signal.argrelextrema, though there are a few different techniques (some on stack overflow) that seem to be just as tricky.

This is most likely related to the comparator function and the waveform.

Long story short:

from scipy.signal import argrelextrema
from matplotlib import pyplot as plt
from scipy import linspace,pi,cos

#Time vector
t = linspace(0,4,4*1000); #4 seconds at Fs=1000Hz
#Phase vector
p = 2.0 * pi * t;
#Signal
S = cos(4.0*p); #A 4Hz sinusoid
# Have a look at it, focus on where the maxima and minima occur.
plt.plot(S);plt.show()
#Let's find some peaks
Pmin = argrelextrema(S, lambda x,y:x<y)
Pmax = argrelextrema(S, lambda x,y:x>y)

Pmin identifies all peaks: [500, 1500, 2499, 3499]

Pmax has skipped some    : [1000,2999]

Try it again, by going for Pmax2 = argrelextrema(S, lambda x,y:x>=y). Pmax2 now is [0,1000,1999,2000,2999,3999]

Now, peak detection is relatively straightforward. The mathematical way to think about it is to get the first derivative and look at where it drops to zero. Every zero of the first derivative is either a local maximum or local minimum. To decide if you are at a peak or a trough, you look at the value of the second derivative at that point. The second derivative of your signal, is the derivative of the first derivative. If that is positive, your candidate point is a local minimum, if it is negative, your candidate point is a local maximum.

That's fine, but costly and there are alternative ways to achieve it. For example, a peak could be defined as the pattern increase, decrease of signal values and a trough could be defined as the pattern decrease, increase of signal values. This is straightforward to understand but in reality, extrema might be "wide". In other words, the pattern of values now is increase, same, same, same, same, same, same, same, decrease. Where is the peak there? (or, conversely, the trough).

You might say, "That's not a problem, the peak is in the middle of the wide top" and that would be a reasonable assumption. In fact, you could achieve this by first "blurring" your data (taking the average) in such a way as to reduce the "wideness" and now you are back in a case of "increase, decrease".

This is achieved in argrelextrema by modifying the order parameter.

So, having said this, why does argrelextrema fails to identify those peaks in the extremely simple case of a cos?

If you take a close look at what the cos is doing, you will observe that:

When the peak is identified:

...,0.99993952,  0.99998889, 0.99999877,  0.99996914...

Did you spot the pattern? (increase,decrease)

When the peak is NOT identified:

..., 0.99995556,  0.99999506, 0.99999506,  0.99995556,...

Did you spot the pattern? (increase,same,decrease)

So, by modifying the comparator function, you enable argrelextrema to widen its search for what an extrema looks like.

Try introducing the <=, >= comparator function.

Hope this helps.


Difference b/w freqz and fft Since both fft ,freqz are used for finding the frequency response of digital filters.
the frequency response using freqz function is given below

and frequency response using fft function is

it seems to be that in case of fft it is mirror image of freqz response.
also i want to find value of freq response at pi/8 how can i find it using any of fft or freqz method.

First, it is instructive to understand the normalised frequency axis of both the fft and the freqz.

By default, in fft, the frequency axis is normalized by the sampling frequency $f_s$. On the other hand, in freqz, the frequency axis is normalised by $f_s/2$, it seems. Having said that, we note that the radian frequency of $\pi/8$ corresponds to the normalized frequency of $\frac{1}{16}=0.0625$ ($2 \pi \frac{f}{f_s}=\frac{\pi}{8}$) in fft and $\frac{1}{8}=0.125$. In  your case, the frequency response at $\pi/8$ seems to be about $-90$ dB.

fft does not provides the mirror image of the freqz. Freqz only shows the half portion of the fft spectrum.


Open access signal processing journal? Signal processing is Elsevier, so it's not open access, and IEEE Transactions on Signal Processing is also not open access. I was wondering if there's any decent journal on par with those that's open access.

I'm trying to submit a paper on discrete time nonlinear filtering.

(there are no tags relating to publication/journal/conference please feel free to add tags or suggest different SE)

The EURASIP Journal on Advances in Signal Processing is a Springer journal with open access.  It sounds like the topic of your paper is appropriate for it.  The impact factor is not as high as IEEE Transactions on Signal Processing though.


How to calculate the mean/center frequency of the spectrum? Assume that the Fourier Transform of $x(t)$ is $X(j\omega)$. And I want to calculate the $\bar{\omega}$ which is the center freqency in the spectrum ( the highest).

And a classmate tells me  this formula $$ \bar{\omega} = \int \omega X(j\omega)\ d\omega$$ could be used to calculate the center frequency ? But it is the first moment, it seems that the result is mean frequency.
Could anyone help me with this question?
What is the mean and center frequency of the spectrum?

There is no solution in Alan's book.

EDIT:
Maybe I did not state the center frequency clearly. I want to find out the frequency where the amplitude spectrum is biggest.
Another situation: assume that the signal is discrete $x[n]$. So we could just compute a period $[-\pi, \pi]$. In this case, this formula could represent the center frequency?

[EDIT] If just want to compute the location of the maximum amplitude, you just have to compute the maximum. You can work in a more generic context, however. As you remarked, the notion of center frequency is not evident. The notion of centroid could be more pertinent.

I assume that your signal is real (since you are talking about a maximum, not well defined in the complex domain). Since its amplitude spectrum is symmetric about $0$, this possible center of mass is useless, since it is the same for all real signals (the zero-frequency). More likely, one is looking at a more useful quantity, based on a one-sided spectrum:

$$\overline{\omega}_W = \int_{\mathbb{R}^+} \omega W(X(\omega))d\omega$$

where $W(X(\omega))$ is a weighting function that (may) depend on the spectrum $X$. The "mean" $\overline{\omega}_W $ will depend on the choice of $W$. A code was given in  the answer to Finding the right measure to compare sound signals in the frequency domain.

You can base it on an absolute value, or something else (positive).  But it is important to have a normalized weight: it should, somehow, sum to one, like for a traditional center of mass. The reason is scaling homogeneity. The mean should not change if you multiply every $X(w)$ by the same quantity.

For orthogonal transforms, based on the energy conservation:

$$ \int [x(t)|^2dt = \int [X(\omega)|^2d\omega $$

it is very common to use a mean/center definition based on an energetic weight, ie
a square of the absolute spectrum. So  every often, the mean frequency is defined as:

$$\overline{\omega}_2 = \frac{\int_{\mathbb{R}^+} \omega|X(\omega)|^2d\omega}{\int_{\mathbb{R}^+}|X(\omega)|^2d\omega}\,.$$

This choice is very sensitive when dealing with time/frequency location uncertainty related to the Weyl-Pauli-Heisenberg inequalities.

But you can use other $L_p$ norms, with:

$$\overline{\omega}_p = \frac{\int_{\mathbb{R}^+} \omega|X(\omega)|^pd\omega}{\int_{\mathbb{R}^+}|X(\omega)|^pd\omega}\,.$$

As $p\to\infty$, this estimator converges to the location of a maximum, as demonstrated in Metric Spaces: Why $L_\infty$ selects the maximum value.

When $p=1$, this is called the spectral centroid, or "brightness" (of sound), or simply "centre of gravity of amplitudes". It is described with other moments in Chapter 6.1 Spectral shape description. When $p=2$, this is sometimes called the power spectral centroid.

I have not found yet useful references for other powers (still searching).


Are delta and double delta features needed when classifying with LSTM I understand that the combination of MFCCs, deltas and double deltas is a good feature to be used with HMMs for keyword detection problems. HMMs are limited by Markov Property and this limitation is overcome by using deltas and double-deltas.

But in the case of LSTMs, they can model long term dependencies very easily. So, doesn't that make the delta and double delta features redundant while using LSTM as the model?

"Very easily" is oversimplification, but overall yes, you do not need to feed deltas into LSTM model.


How does cascading low pass modules affect pass band ripples? This question 10.21 is taken from GATE IN paper 2015.

I don't understand the answer given,

I mean why is the frequency response at the pass band frequency considered to be $(1-\delta)$, isn't the ripple equally distributed on both sides, so shouldn't the average be One?

Finally why is effective ripple calculated by subtracting from one?

The passband region of an equi-ripple filter is the region where its magnitude is in the interval $[1-\delta,1+\delta]$, where $\delta$ is the maximum approximation error. At the passband edge the value of the magnitude must equal $1-\delta$ (as shown in your drawing). The passband edge is the last frequency where the magnitude is still inside the interval $[1-\delta,1+\delta]$, and since the magnitude then falls off towards the stopband, the magnitude at the passband edge must be at the lower limit of the passband magnitude, which equals $1-\delta$.

The total magnitude response of a cascade of $M$ identical filters is

$$|H_{total}(\omega)|=|H(\omega)|^M$$

Consequently, the total magnitude at the passband edge is $(1-\delta)^M$. This value is then equal to $1-\delta_{total}$, which gives

$$\delta_{total}=1-(1-\delta)^M$$


What is meant by *sampling* in terms of the *sampling theorem*? Let $y:\left(-\frac T2,\frac T2\right)\to\mathbb{C}$ be a square integrable function. The Fourier coefficients of $y$ are $$\underline{Y}(k):=\frac 1T\int_{-T/2}^{T/2}y(t)e^{-i\omega_kt}\;dt\;\;\;\text{with }\omega_k:=k\frac{2\pi}T$$
for $k\in\mathbb{Z}$. The Fourier polynom of degree $n\in\mathbb{N}$ of $y$ is $$\mathcal{F}^{-1}_n[y](t):=\sum_{k=-n}^n\underline{Y}(k)e^{i\omega_kt}$$ and $$\mathcal{F}^{-1}[y]:=\lim_{n\to\infty}\mathcal{F}_n^{-1}[x]$$ is called inverse Fourier transformation of $y$. Now, I've got two questions:

What is meant by sampling (in terms of the sampling theorem)? From my understanding, if we know the period $T$ all we need to "store" are the values $\underline{Y}(k)$. We cannot store all values, so we need to choose a "huge enough" $n$ and store only the values $\underline{Y}(-n),\cdots,\underline{Y}(n)$. So, where does "sampling" come into play? The only thing I could imagine is numerical integration: We consider an equidistant grid $$x_j=\left(\frac jN-1\right)\frac T2\;\;\;\text{for }j=0,\ldots,2N$$ and apprximate $\underline{Y}(k)$ using the composite trapezoidal rule, i.e. $$\underline{Y}(k)\approx\frac{1}{2N}\sum_{j=0}^{2N-1}y\left(x_j\right)e^{-i\omega_kx_j}$$ By doing so, we didn't take the whole "signal" $y$, but only the "sample points"$\left(x_j,y\left(x_j\right)\right)$ into account. Is this meant by "sampling"?
Does the sampling theorem make a statement about $n$ or $N$ or something else?

I think, with respect to sampled function, that i disagree a bit with MBaz.  perhaps i misunderstand what he/she says regarding needing an infinite number of coefficients or samples for this special case of sampled periodic functions.  we had a little discussion about this at comp.dsp, and i'm gonna make use of $\LaTeX$ to spell out the same points.

of course, whether $x(t)$ is periodic or not, if $x(t)$ is real and is sufficiently bandlimited (there are no frequency components as high or higher than $\frac{f_\text{s}}{2} = \frac{1}{2T_\text{s}}$), then the samples:

$$ x[n] \triangleq x(n T_\text{s}) $$

are sufficient to completely represent the continuous-time $x(t)$.  if $x(t)$ never repeats, than an infinite number of discrete $x[n]$ are necessary to represent $x(t)$ for all $t$.

but if $x(t)$ is periodic,

$$ x[n+N]=x[n] \quad\quad \forall n,N \in \mathbb{Z} $$
and
$$ x(t+T)=x(t) \quad\quad \forall t $$

then samples existing over the span of one period are sufficient to represent $x[n]$ and $x(t)$.  in order for the $x[n]$ to be periodic in the sampled domain:

$$ \begin{align}
x[n+N] & = x\left((n+N)T_\text{s} \right) \\
& = x\left(nT_\text{s}+NT_\text{s} \right) \\
& = x\left(nT_\text{s}+T \right) \\
\end{align} $$
because $ x(t+T) = x(t) $.

which means the obvious, your function period $T$ has to be the same as $N$ times your sampling period $T_\text{s}$.

$$ T = N T_\text{s} $$

now, what we know about this sampled periodic function is that $N$ samples of $x[n]$ are sufficient to tell us all about $x[n]$, and since $x(t)$ is bandlimited, $x[n]$ and the $N$ samples that fully define it, are sufficient to fully describe $x(t)$.

if $x(t)$ is sufficiently bandlimited (as above), then

$$ \begin{align}
x(t) & = \sum\limits_{n=\infty}^{+\infty} x[n] \operatorname{sinc}\left(\frac{t-nT_\text{s}}{T_\text{s}}\right) \\
& = \sum\limits_{n=\infty}^{+\infty} x[n] \frac{\sin\left(\pi\frac{t-nT_\text{s}}{T_\text{s}}\right)}{\pi\frac{t-nT_\text{s}}{T_\text{s}}} \\
& = \sum\limits_{m=\infty}^{+\infty} \sum\limits_{n=0}^{N-1} x[n+mN] \frac{\sin\left(\pi\frac{t-(n+mN)T_\text{s}}{T_\text{s}}\right)}{\pi\frac{t-(n+mN)T_\text{s}}{T_\text{s}}} \\
& = \sum\limits_{n=0}^{N-1} \sum\limits_{m=\infty}^{+\infty} x[n] \frac{(-1)^{mN}\sin\left(\pi\frac{t-nT_\text{s}}{T_\text{s}}\right)}{\pi\frac{t-(n+mN)T_\text{s}}{T_\text{s}}} \\
& = \sum\limits_{n=0}^{N-1} x[n] \sin\left(\pi\frac{t-nT_\text{s}}{T_\text{s}}\right) \sum\limits_{m=\infty}^{+\infty} \frac{\frac{T_\text{s}}{\pi}(-1)^{mN}}{t-(n+mN)T_\text{s}} \\
\end{align} $$

$$ x(t) =
\begin{cases}
\sum\limits_{n=0}^{N-1} x[n] \frac{\sin\left(\pi\frac{t-nT_\text{s}}{T_\text{s}}\right)}{N\tan\left(\frac{\pi}{N}\frac{t-nT_\text{s}}{T_\text{s}}\right)}, & \text{if }N\text{ is even} \\
\sum\limits_{n=0}^{N-1} x[n] \frac{\sin\left(\pi\frac{t-nT_\text{s}}{T_\text{s}}\right)}{N\sin\left(\frac{\pi}{N}\frac{t-nT_\text{s}}{T_\text{s}}\right)}, & \text{if }N\text{ is odd}
\end{cases} $$

proving the latter takes a little bit.  you might recognize the $N$ odd case as the Dirichlet_kernel.  the $N$ even case looks a teeny bit different.  but both show exactly how the $N$ samples that fully define the sampled and periodic $x(t)$ are combined to get $x(t)$.

now, since $x(t)$ is also periodic with period $T$, then

$$ \begin{align}
x(t) & = x\left(t+T \right) \\
& = x\left(t+N T_\text{s} \right) \\
& = \sum\limits_{k=-\infty}^{+\infty} X[k] e^{i 2 \pi \frac{k}{T} t} \\
& = \sum\limits_{k=-\infty}^{+\infty} X[k] e^{i 2 \pi \frac{k}{N T_\text{s}} t} \\
& = \sum\limits_{k=-\lfloor \frac{N}{2} \rfloor}^{+\lfloor \frac{N}{2} \rfloor} X[k] e^{i 2 \pi \frac{k}{N} \frac{t}{T_\text{s}}}
\end{align} $$

where $\lfloor\cdot\rfloor$ is the $\operatorname{floor}(\cdot)$ operator and

$$ \begin{align}
X[k] & = \frac{1}{N} \sum\limits_{n=0}^{N-1} x[n] e^{-i 2 \pi \frac{nk}{N}} \\
& = \mathcal{DFT}\{x[n]\}
\end{align} $$

there's actually something to fudge (a factor of $\frac{1}{2}$) about $X\left[\frac{N}{2}\right]$ for the $N$ even case:

$$ \begin{align}
X\left[-\frac{N}{2}\right] = X\left[\frac{N}{2}\right] & = \frac{1}{2} \ \frac{1}{N} \sum\limits_{n=0}^{N-1} x[n] e^{-i 2 \pi n\frac{n(N/2)}{N}} \\
& = \frac{1}{2N} \sum\limits_{n=0}^{N-1} x[n] (-1)^n \\
\end{align} $$

note that:

$$ \begin{align}
x[n] = x(n T_\text{s}) & = \sum\limits_{k=-\lfloor \frac{N}{2} \rfloor}^{+\lfloor \frac{N}{2} \rfloor} X[k] e^{i 2 \pi \frac{k}{N} \frac{t}{T_\text{s}}}\bigg|_{t=n T_\text{s}} \\
& = \sum\limits_{k=-\lfloor \frac{N}{2} \rfloor}^{+\lfloor \frac{N}{2} \rfloor} X[k] e^{i 2 \pi \frac{k}{N}n} \\
& = \sum\limits_{k=0}^{N-1} X[k] e^{i 2 \pi \frac{nk}{N}} \\
& = \mathcal{iDFT}\{X[k]\}
\end{align} $$

if you deal with that $\frac{1}{2}$ fudging for $X\left[\frac{N}{2}\right]$ for the $N$ even case.  this is because $X[k+N]=X[k]$ for all $k$.

a periodic continuous-time function can be described with a countable set of Fourier coefficients.  a bandlimited periodic continuous function can be described with a finite set of $N$ Fourier coefficients, just as well as it can be described with a finite set of $N$ samples.


How to cut off the initial disturbance automatically? I would like to find a way to cut off the initial disturbance (e.g. the part tf_n<20). Though it is easy to do manually for single figure, I need a way to automatically identify the steady point, as I need to do it 100 times otherwise.

I would prefer to do it by Matlab.

I just come up with a solution when I am writing question...

I can compare the amplitudes of this period and last steady periods (for my case there is a fixed responding frequency). When the relative error is small enough for 5 or 10 times, it can be considered as a good cutoff point.

below a Matlab script I used, and it worked well for my case.

function [t1,y1,t2,y2]=autocut(f1,t1,y1,t2,y2)
    errlim=10e-2;
    cut1=findcut(f1,t1,y1,errlim);
    cut2=findcut(f1,t2,y2,errlim);
    cutp=min([cut1 cut2]);
    t1=t1(cutp:end);
    t2=t2(cutp:end);
    y1=y1(cutp:end);
    y2=y2(cutp:end);
end

function cutnt=findcut(f1,t,y,errlim)
dt=t(20)-t(19);
T1n=ceil(1/f1/dt)*2;%T1n points in one period, number of periods for A=(max-min)/2
tn=length(t);
nperiod=floor(tn/(T1n));
%% Ast is the amplitude of last 5 periods with steady vibration
tstartn=(nperiod-5-1)*T1n+1;tendn=nperiod*T1n;
Ast=(max(y(tstartn:tendn))-min(y(tstartn:tendn)))*0.5;
%% if the relative error is continuously small for 15 periods,
s=0;%count periods that satisfies the condition, 5 statisfaction periods define the steady
for i1=1:nperiod %
    tstartn=(i1-1)*T1n+1;
    tendn=i1*T1n;
    yamp(i1)=(max(y(tstartn:tendn))-min(y(tstartn:tendn)))*0.5;
    err(i1)=abs((yamp(i1)-Ast)/yamp(i1));%relative error
    if err(i1)<errlim
        s=s+1;
        if s==5
            break
        end
    end
end

cutnt=tendn;


Determining which Filter from a Z-Plane Plots? How do i determine which FIR filter (LP, HP, BS, BP) it is from looking at it's z-plane plot?

just to add a little to Matt's answer.  if you know where your poles, $p_n$, and zeros, $q_m$, are your transfer function is

$$ H(z) = H_0 \frac{(z-q_1)(z-q_2)(z-q_3)...(z-q_M)}{(z-p_1)(z-p_2)(z-p_3)...(z-p_N)} \quad \quad M \le N $$

the complex frequency response (which contains both magnitude and phase information) is

$$ H(e^{j\omega}) = H_0 \frac{(e^{j\omega}-q_1)(e^{j\omega}-q_2)(e^{j\omega}-q_3)...(e^{j\omega}-q_M)}{(e^{j\omega}-p_1)(e^{j\omega}-p_2)(e^{j\omega}-p_3)...(e^{j\omega}-p_N)} $$

and the magnitude response (which is where you determine if it's a LPF or BPF or HPF or something else like that) is:

$$ |H(e^{j\omega})| = |H_0| \frac{|e^{j\omega}-q_1| \ |e^{j\omega}-q_2| \ |e^{j\omega}-q_3|...|e^{j\omega}-q_M|}{|e^{j\omega}-p_1| \ |e^{j\omega}-p_2| \ |e^{j\omega}-p_3|...|e^{j\omega}-p_N|} $$

and $e^{j\omega}$ is a point on the unit circle.  what this means is, as your frequency increases from DC (or $\omega=0$), that corresponds to moving along the unit circle from $z=e^{j0}=1$ to other points on the unit circle.  then as your point on the unit circle gets close to a zero $q_m$, a factor of $|e^{j\omega}-q_m|$ gets small.  $|e^{j\omega}-q_m|$ is literally the distance your point on the unit circle is from that particular zero.

likewise for the poles, but your poles are in the denominator, so as your point on the unit circle gets close to a pole $p_n$, a factor of $|e^{j\omega}-p_n|$ gets small, but since you're dividing by that small number, the effect on the gain magnitude is to increase it.

so a LPF will have poles closer to the DC point of $z=e^{j0}=1$ and, perhaps, zeros closer to Nyquist $z=e^{j\pi}=-1$ so when you're close to DC you're multiplying by big numbers and when you're close to Nyquist, you're multiplying by small numbers.

likewise for a HPF, you'll have poles closer to Nyquist and zeros closer to DC.

for a BPF, you might have zeros close to both DC and Nyquist and poles close to the unit circle at other places where they will boost the gain for those bandpass frequencies.

that's how you look at the pole-zero plots on the $z$-plane to judge if it's a LPF or whatever.


Low PAPR performance for Monochromatic frequency signals I understand what is Monochromatic frequency signals from here, but could somebody kindly give a brief explanation about "why does a monochromatic frequency signal result in poor PAPR (Peak to average Power Ratio) performance"?.
(Note- topic for achieving Good PAPR(constant Envelope) for Radars)

The term monochromatic signal is really just a more complicated way of saying a sinusoid, a function that only contains spectral content at one frequency. Peak-to-average power ratio (PAPR) is a bit of a mushy term if you start to look for a theoretical definition, as power by definition is the time rate of delivery of energy over some period of time; the measurement period that you choose will affect the power value that you select.

For our purposes, let's just assume that the "instantaneous" delivered power is proportional to the square of the signal's amplitude (as if it was driving a resistive load). Let's assume the peak power of the signal, then, corresponds with the peak value of the sinusoidal signal $x(t)$:

$$
x(t) = A \cos{\left(2\pi f t \right)}
$$

$$
P_{peak} \propto A^2
$$

To find the average power delivered by the sinusoid, we use the RMS value of a sinusoidal signal:

$$
A_{rms} = \frac{A}{\sqrt{2}}
$$

$$
P_{av} \propto \frac{A^2}{2}
$$

The coefficient of proportionality is the same for both cases, so the ratio is:

$$
\text{PAPR} = \frac{A^2}{\frac{A^2}{2}} = 2
$$

This PAPR actually isn't that bad when compared to many modulated signal types. As an example, OFDM signals often have PAPRs of 10 dB or more, which can complicate the design of amplifiers that are suitable for driving said signals for transmission.


Non-decaying Harmonics in Vibration Signal I am analyzing vibration data collected by an accelerometer positioned in a stamp press machine. The fundamental frequency is around 10 Hz. The problem here is that when I take the FFT of the data I see harmonics of the fundamental frequency until 5000 Hz every 10 Hz (of course). I would expect some visible harmonics in low-frequency but with a decaying behavior until, at some point, it disappears. I attached 2 figures. The first one with a zoom from 0 to 1000 Hz and the other from 4000-5000 hz. Note that the Harmonics dominate the signal.

And here is a sample of the data:

I have being trying to convince myself that is normal but I can´t. What do you guys think? Can it be an effect of something that I should discard?

Thank you!

Well, looking at your time-domain signal, it looks a bit similar to a pulse train. And the Fourier Transform of a Pulse train with distance T is a Pulse train with distance 1/T

Hence, if you see a pulse every 10Hz, this corresponds to pulses every 0.1seconds in the time domain. That's what the time-domain signal looks like.


Does sampling in the frequency domain cause time-domain aliasing? Let's say I have an impulse response $h[n]$.

I analyze the power spectrum of that impulse response similar to fourier transformed $h[n]$ corresponding to roughly $H[f]$.

Now I compare $H[f]$ with some target $H_{0}[f]$ by subtracting $H_{0}[f]-H[f]$.

I create a compensation filter based on this via inverse fourier transform: $\mathrm{ifft}(H_{0}[f]-H[f])$.

But this is inadequate because time-domain aliasing is observed: Why?

What I want to understand is, just as sampling in the time-domain cause aliasing in the frequency domain, is this operation of creating a filter via ifft akin to sampling in the frequency domain and likewise causing aliasing in the time domain?

Or is this because the $H[f]$ was not sufficiently padded and the result is a circular convolution that causes ringing (but is this time-domain aliasing?)

Is this a direct result of the duality of the DFT? In that, just as a sinusoid of a given frequency will show up as a frequency-axis shifted impulse in the frequency domain, an envelope of some frequency will show up as a time-axis shifted impulse in the time domain?

Confusedly yours,
@panthyon

The answer is: yes, sampling in the frequency domain causes aliasing in the time domain, exactly like the dual case: sampling in the time domain causes aliasing in the frequency domain.

There are many ways to see this. One standard way is to sample the discrete-time Fourier transform (DTFT) of a discrete-time signal by multiplying it with a Dirac comb and see what this does to the time domain signal. Here, however, I would like to show it in another way, namely by taking samples of the DTFT and asking which periodic signal has these samples as its discrete Fourier transform (DFT).

Let $Y(e^{j\omega})$ be the DTFT of the sequence $y[n]$ (note that $y[n]$ is not necessarily of finite length):

$$Y(e^{j\omega})=\sum_{n=-\infty}^{\infty}y[n]e^{-jn\omega}\tag{1}$$

Now we take samples of $Y(e^{j\omega})$ at $\omega_k=2\pi k/N$, $k=0,1,\ldots,N-1$, and derive the time domain sequence with a DFT equal to these samples:

$$\begin{align}Y_k\doteq Y(e^{j2\pi k/N})&=\sum_{n=-\infty}^{\infty}y[n]e^{-j2\pi nk/N}\\&=\sum_{l=-\infty}^{\infty}\sum_{n=0}^{N-1}y[n+lN]e^{-j2\pi nk/N}\\&=
\sum_{n=0}^{N-1}\sum_{l=-\infty}^{\infty}y[n+lN]e^{-j2\pi nk/N}\\&=\sum_{n=0}^{N-1}\tilde{y}[n]e^{-j2\pi nk/N}\tag{2}\end{align}$$

where $\tilde{y}[n]$ is an aliased version of the original sequence $y[n]$:

$$\tilde{y}[n]=\sum_{l=-\infty}^{\infty}y[n+lN]\tag{3}$$

Note that the right-hand side of $(2)$ is the DFT of $\tilde{y}[n]$.
Eq. $(2)$ shows that the samples $Y_k$ of the DTFT $Y(e^{j\omega})$ are the DFT coefficients of $\tilde{y}[n]$, which is an aliased version of $y[n]$. Consequently, sampling in the frequency domain causes time-domain aliasing, which is obviously the dual of the well-known aliasing in the frequency domain caused by sampling in the time domain.


Is there a noise cancellation technique for over-lapping spectra? I understand Wiener filter is used to reduce noise, Gaussian noise as well, once we give the noise characteristic. One may derive the filter coefficient in adaptive way as well But can it remove noise if the noise is in the same band as the desired signal(EDIT: without loosin in-band energy of the signal)? Is there an such filter which is applicable for in band noise?  As per my understanding there is none and is not possible.
Edit: The signal is not assumed to be stationary.

Any suggestion, comment please?

First of all the Wiener filter does not remove the noise but reduce it for WSS signals.

It does this based on the relationship between the power spectral density of the clean signal $x[n]$ and imposed noise $v[n]$, both assumed as WSS. The frequency response of a noncausal IIR Wiener filter is given as

$$
H(\omega) = \frac{ P_x(\omega) }{ P_x(\omega) + P_v(\omega) }
$$

where $P_x$ and $P_v$ are PSDs of clean signal and uncorelated noise imposed on it. As can be seen, the filter spectrum involves all the bands including where noise and signal spectrums are overlapping too, and that's one of the benefits of using such an advanced tool.

Note that this is not a non-destructive noise reduction and it wont necessarily preserve the waveform of the clean signal, but its PSD is improved; made closer to clean signal's PSD.