Using the given identities, find the inverse DTFT Using the given identities,

$$ a^nu[n] \Longleftrightarrow \frac{1}{(1-ae^{-jw})}$$
and
$$\delta[n-k]\Longleftrightarrow e^{-jwk}$$

Find the inverse DTFT of,

$$ H(e^{jw}) = B\cdot\frac{e^{-jw}}{(1-ae^{-jw})}$$

my attempt:

$$h[n] = B\cdot\delta[n-1]a^nu[n]$$

It seems straight forward enough, just plug in the inverse dtft. This is not correct though, there is no delta to be found in the correct solution.

The correct answer is:

$$h[n] = Ba^{n-1}u[n-1]$$

It's like the delta disappeared and the delta convolution property was used.

$$x[n]*\delta[n-1] = x[n-1]$$

But using the identities above, I do not want to use convolution. I am confused why the delta disappeared and the signals it was multiplied with became shifted, any help in understanding how to get the correct answer using the identities would be appreciated!

As pointed out in the comments, it is important for you to know and understand that multiplication in one domain (time or frequency) corresponds to convolution in the other domain. And, as a consequence, multiplication with $e^{-j\omega}$ in the Fourier domain corresponds to a delay of one sample in the time domain (i.e., convolution with $\delta[n-1]$). So you have

$$\begin{align}\mathcal{F}^{-1}\left\{B\frac{e^{-j\omega}}{1-ae^{-j\omega}}\right\}[n]&=B\cdot\mathcal{F}^{-1}\left\{\frac{e^{-j\omega}}{1-ae^{-j\omega}}\right\}[n]\quad (\textrm{due to linearity of }\mathcal{F}^{-1})\\&=B\cdot\mathcal{F}^{-1}\left\{\frac{1}{1-ae^{-j\omega}}\right\}[n-1]\\&=B\cdot a^{n-1}u[n-1]\end{align}$$

where I used $\mathcal{F}^{-1}$ to denote the inverse discrete-time Fourier transform.


Approximate the change in perceived gain from applying transfer function to audio signal I have an audio signal that I am applying an analogue transfer function $H(s)$ to. I would like to approximate the change in gain ($\Delta \text{ Gain}$) as perceived by the listener for an arbitrary audio signal sent through the filter $H(s)$.

I was thinking something along the lines of:

$$
\Delta \text{Gain} = \frac {\displaystyle\int_{0}^{\pi} \lvert H(\omega)\rvert d\omega} {\pi}
$$

This is based on an assumption that the formula above gives a sensible approximation of gain across the spectrum by calculating the averaging gain.

Is my thinking sensible or is there some proper theory that will give a better answer that I can refer to?

You want to calculate loudness difference between input and output.

Loudness depends on the power spectral density of a signal; knowing only the transfer function of your filter is not enough for calculating PSD.

For example, suppose your filter is a lowpass filter, with attenuation 20 dB in the stop band. If your input signal is completely contained below the cutoff frequency, then no perceived gain (input and output are similar), whereas if your input signal is all above the cutoff frequency, then you will perceive a -20dB gain.

Related question on calculating loudness.


SNR calculation confusion Follow up to question and comments posted here

The SNR calculation method posted as answers were implemented in matlab

Implementation of Deve's answer:

% the starting and ending point of speech in samples obtain through Audacity;
% select a portion of the waveform where speech is seen to be present
begin_speech = 25400; end_speech = 26800;

snr_before = mean( data_s(begin_speech:end_speech) .^2) / mean( data_n(begin_speech:end_speech) .^2);
db_snr_before = 10*log10( snr_before );     % same thing, but in dB

% calculation on data after noise reduction follows:
residual_noise = data_s(begin_speech:end_speech) - output(begin_speech:end_speech);
snr_after = mean( data_s(begin_speech:end_speech) .^ 2)/mean( residual_noise .^ 2);
db_snr_after = 10*log10( snr_after );

The result is

db_snr_before = 15.6853
db_snr_after = -22.0388

Implementation of dspGuru's answer:

begin_speech = 25400; end_speech = 26800;

% before noise reduction
DB1 = 10*log10( var(data_n(begin_speech:end_speech)));
DB2 = 10*log10( var(data_s(begin_speech:end_speech)));
db_SNR_before = DB2 - DB1;

% after noise reduction
[data_stnr, fs_stnr, nbits_stnr] = wavread('AdMStry_thru_nr.wav');
Z = output(begin_speech:end_speech) - data_stnr(begin_speech:end_speech);
db_residual_noise = 10*log10( var(Z) );
db_speech = 10*log10( var(data_stnr(begin_speech:end_speech)));
db_SNR_after = db_speech - db_residual_noise;

The result is :

db_SNR_before = 15.6848
db_SNR_after = 0.9252

The begin_speech and end_speech represent those parts of the wavfile, in number of samples, where there is presence of speech for certain (obtained through audacity by observing the waveform directly)

Why is the SNR less after passing through noise reduction?

In the second method (dspGuru), to calculate SNR after noise reduction, to get the Signal value on the numerator, I have passed just the pure speech through the noise reduction algorithm as well. If I do the same thing in the first (Deve) method, the results of the two methods become almost equal

Why is this happening? Why is SNR after less than SNR_before ?

Here is a screenshot of the waveforms:

As Jason R has already mentioned in his comment, the noise reduction algorithm obviously introduces some amplification which renders the SNR estimates we proposed in our earlier answers useless. This becomes clear if we look at how the SNR is calculated.

Let $x_k$ be the speech signal without noise and let $n_k$ be the noise. Then the SNR before noise reduction is
$$
\gamma_{before} = \frac{\sum x_k^2}{\sum n_k^2}
$$

Now the noise reduction (NR) is applied to both signals yielding the pure speech signal after NR $\tilde{y}_k \neq x_k$.  Applying NR to the noisy signal produces some new signal $y_k + \tilde{n}_k$. As can be seen from your screenshots $y_k \neq \tilde{y}_k$. If we apply dspGurus method for calculating the SNR after NR we obtain
$$
\gamma_{after} = \frac{\sum \tilde{y}_k}{\sum [y_k - \tilde{y}_k + \tilde n_k]^2}
$$
The denominator not only conatins the residual noise but also an additional error $y_k - \tilde{y}_k$. This leads to a wrong estimation and a much lower SNR value.

To overcome this problem, you could try to normalize the signals. This will be difficult, because obiously the NR algorithm scales the noisy and the noise free signal differently. If you've implemented the algorithm yourself, you should try to fix this.

Another possible solution is to estimate the noise power in a region where the original signal contained silence and to calculate the total signal power in a region that contains noise and signal:

Take a region from the noisy signal after NR that only contains noise and calculate the mean power denoted as $\sigma_{n}^2$
Take a region from the noisy signal after NR that contains noise and speech and calculate the mean power denoted as $\sigma_{n,s}^2$

Now the SNR after NR can be estimated as
$$
\gamma_{after} = \frac{\sigma_{n,s}^2 - \sigma_{n}^2}{\sigma_{n}^2}
$$

How good this estimation is, depends on how good the following prerequisites are fulfilled in your case:

A large number of samples should be used to calculate the mean powers. You used 1401, that's not much.
Noise and signal after NR have to be uncorrelated. The NR algorithm possibly introduces noise that is a function of the signal. So this condition might not be fulfilled.
The mean noise power should be constant over time. I think this is fulfilled.


Types of interpolation used for reconstruction in DSP? What are the different types of interpolation used in DSP for reconstruction of analog signal from discrete/digital signal

I am able to somehow learn two types of interpolation
1st is "zero order hold interpolation"and 2nd is "linear interpolation(first order interpolation)"

Reconstructed output looks like stair case or square wave for zero order hold interpolation while reconstructed output for linear interpolation(1st order interpolation) looks like a triangular wave

There are other types of interpolation besides above?

What are those types?  (That types give a much better/smoother reconstructed signal as compared to two types i discussed above)

What is the order of polynomial in those types (just like "zero order hold" type had zero order and "linear interpolation" type had order one)

For those other types of interpolation, what/how are the appearances/looks of reconstructed output waves ? (Please kindly try to include those wave graphs in your answer)

By reconstruction, i mean conversion of digital into analog signal

Zero-order hold will result in a piecewise-constant waveform.

Linear interpolation will result in a piecewise-linear waveform.

If you want a piecewise-quadratic or piecewise-cubic or higher order polynomial interpolation, it will not appear much different from the original bandlimited waveform.


An equality in DOA estimation Suppose we use a uniform linear array for DOA estimation. Assuming neither of the signals and noise are correlated, and the noise power on each sensor is equal, it's well known that the covariance of the received signal can be decomposed as:
$$\boldsymbol{R}=\boldsymbol{A}\boldsymbol{P}\boldsymbol{A}^*+\sigma^2 \boldsymbol{I}=\boldsymbol{S}\boldsymbol{\varLambda_{\mathrm{s}}}\boldsymbol{S}^*+\sigma^2\boldsymbol{G}\boldsymbol{G}^*,$$
where the eigenvectors in $\boldsymbol{S}$ and $\boldsymbol{G}$ are selected to be orthogonal and normalized.
My question is how to prove the equality $$\boldsymbol{G}\boldsymbol{G}^*=\boldsymbol{I}-\boldsymbol{A}(\boldsymbol{A}^*\boldsymbol{A})^{-1}\boldsymbol{A}^*$$ in the paper Maximum Likelihood Methods for Direction-of-Arrival Estimation. It seems very basic but I can't come up with the proof.

I think I know the answer now. If I made an error, please correct me.

Assume there are $m$ sensors and $d$ sources. As the eigenvectors in $S$ and $G$ of the Hermitian matrix $R$ is selected to be orthogonal and normalized,
$$UU^*=\begin{bmatrix}S&G\end{bmatrix}\begin{bmatrix}S^*\\G^*\end{bmatrix}=SS^*+GG^*=UU^{-1}=I_m;$$
$$S^*S=I_d;$$
$$\mathrm{span}\{G\}=\mathrm{span}\{S\}^{\bot}.$$

We know the noise subspace $\mathrm{span}\{G\}$ and $\mathrm{span}\{A\}$ are orthogonal as in the MUSIC algorithm:
$$\mathrm{span}\{G\}=\mathrm{span}\{A\}^{\bot}.$$
So $\mathrm{span}\{A\}=\mathrm{span}\{S\}$, which implies their projection matrices into this subspace are equal:
$$A(A^*A)^{-1}A^*=S(S^*S)^{-1}S^*.$$
Using $S^*S=I_d$, we get
$$A(A^*A)^{-1}A^*=SS^*;$$
with $SS^*+GG^*=I_m$, we get
$$GG^*=I_m-SS^*=I_m-A(A^*A)^{-1}A^*.$$


Use GNU Radio to extract information from a VOR signal I want to design a GNU Radio project to extract the information of a VOR signal, then find the angle to the north and the Morse code of the station.

Is there any example that shows how to do it?

I'm not aware of any existing application.
Azimuth Estimation
However, wouldn't VOR (VHF Ominidirectional Ranging) consist of two relevant signals:

the reference signal, a 30 Hz signal frequency-modulated onto a  carrier not too far of the nominal system frequency
the rotating beam signal, which is phase-delayed proportionally to the azimuth

with the direction information in the phase difference between reference signal and beam signal?
so your first implementation would look something like

demodulate the FM signal and estimate the phase
estimate the phase of the pulsed beam signal at its maximum
compare these two

A more stable second implementation would phase-lock an internal oscillator to the reference signal. An even more advanced receiver could then e.g. apply a Kalman filter to the observed azimuth.
Morse Decoding
I'm also not aware of any Morse/CW receiver in existence (there probably is – a good idea is just asking for this isolatedly on http://ham.stackexchange.com , or the GNU Radio mailing list.
Decoding human transmissions is probably non-trivial, since humans tend to have suboptimally stable symbol clocks, but an automatic transponder might be relatively easy to decode – just find the shortest reasonable "dit", everything significantly longer is a "dah".


Minimum phase FIR method I am trying to make a minimum phase filter (in wxMaxima) according to these steps:

first create a "normal" FIR (a simple sinc, wc=0.4, random example, but in the pictures a remez with Octave)
calculate the zeroes (allroots() function)
delete the zeroes outside the unit circle (anything>1/0.98, tolerances&co)
delete half of the zeroes on the unit circle (the conjugate pairs, 0.98 < x < 1/0.98)
reconstruct the polynomial => coefficients.

As I understand it, these are the correct steps to follow, but this is what I'm getting:

and this only if I don't delete half the zeroes. If I do, no matter how I sort, choose, etc, I get bogus results. Is this the right way to achieve a minimum phase? I can provide the wxMaxima approach if needed (I tested it first to see if I can reconstruct the original impulse and I could do it, so obtaining the coefficients is not the problem, but the how to).

I wanted to add this later that day but got busy. Here's a graph showing the frequency response of the four filters, ignore the blatant uglyness in attenuation, it's meant to be an example. The 2nd is by zero inversion (and doubling), the 3rd is by spectral factorization and the 4th is like the 3rd only done by geometric mean between adjacent zeroes, to avoid the manual tweaking of adding the DC offset to h[M/2]:

Here are the zeroes:

The idea with the geometric mean came after realizing that the response of the equiripple wasn't quite equiripple, so adding $\delta_2$ to h[M/2] was tricky (you can see the green zeroes aren't quite overlapping). This also gave me the idea that I can use it with non-equiripple, too (that pdf implies the method can only be applied there). Sure enough, here's the result:

Of course, windowing with a Kaiser (Dolph-Chebyshev, or exponential, or cosh) window(s) should be necessary to compensate for the stop-band. Here's a quick example with a Hamming window:

I think the geometric mean ($\sqrt{z_1 \cdot z_2}$) comes as a good alternative, as long as the zeroes are fairly evenly-spaced, otherwise choosing between non-overlapping zeroes (see the green ones) results in going away from their center.

After this I reached the conclusion that if you need a minimum-phase FIR then it's safe to go for an inverse Chebyshev IIR for the sake of less computation power. The half-order minimum phase might be worth it in terms of delay, but it's a pain matching the amplitude. Here are some results with LTspice, showing two FIRs with N=64, linear and nimimum phase, and an inverse Chebyshev IIR (not Cauer, chosen for flatness and similarity in response) with N=11 and matching attenuation:

And, for the sake of comparison, here are the phases of the min.-ph. FIR and the inv. Cheb. IIR:

...and the step responses of all three:

You made a minimum-phase filter but with a different magnitude response than the original linear phase filter. What you have to do to keep the magnitude the same is to reflect the zeros outside the unit circle into the circle. Since for linear phase filters the zeros are symmetrical with respect to the unit circle (or on the circle), you simply have to double each zero inside the unit circle and leave the zeros on the circle unchanged. Apart from doubling each zero inside the circle you also need to apply a scale factor $1/z_i$ for each zero that gets reflected (where $z_i$ is the location of the zero).

So each pair of zeros of the linear phase filter (assuming $|z_i|<1$)

$$(1-z_iz^{-1})(1-z^{-1}/z_i^*)$$

needs to be replaced by the following scaled double zero:

$$\frac{1}{|z_i|}(1-z_iz^{-1})^2$$

This will keep the magnitude unchanged while transforming the linear phase function to a minimum-phase function.


Choosing the length of a signal when calculating the FFT I'm trying to transforming an already know signal from time domain to frequency domain using fft function in MATLAB. And I got different values of amplitude when I change the value of length of signal $L$. I'm therefore wondering if there is a rule to choose the length of signal when calculating FFT.

I have an acoustic signal got from a sensor and need to analyze it (I have no idea about its frequency range, amplitude, and so on.) I want to transform it into frequency domain using fft function in MATLAB. And I don't know how can I choose the length of signal.

Could somebody please help?

Mai,

The length of the FFT depends on what application you are doing. A very course summary follows:

Same size FFT:

Analysis: This just means you want to 'analyse' the signal - look and see what type of spectrum it has, maybe patterns in the spectrums, etc. The 'usual' default is to simply use the FFT length equal to the length of your signal. Example: "What is the frequency content of my signal look like?"z = abs(fft(x));
Circular convolution: In circular convolution, you want to convolve two signals circularly, as opposed to 'normal' linear convolution. So if you have a signal $x$ of length $L$, and you want to circularly convolve with signal $y$ of length $L$ also, you can multiply their same length FFTs with each other, and then take the IFFT. Example: "What is the circular convolution of $x$ and $y$?"
z = ifft(fft(x).*fft(y))

Different Length FFT:

Analysis: You can also use different (bigger) size FFTs for analysis as well. For example, if your signal $x$ was of length $L$, you can still analyze its frequency content choosing an FFT size of length, say, $10*L$. If you do this, you are actually interpolating the frequency spectrum. This just means that the spectrum looks more 'smooth', and in some cases you are able to better see where peaks might exist. Importantly though, this is not the same as increasing your frequency resolution. If someone tells you to "increase frequency resolution by taking a bigger size FFT", tell them they are wrong, and that this just interpolates the frequency domain. (Unless of course that person is your boss in which case you have to solve the non-linear optimization problem of telling them they are wrong without actually telling them they are wrong :-) ). (To actually increase frequency resolution, you need to have more samples of your original time domain signal). Anyway, the example motivation for this might be: "We analyzed our time domain signal using same size FFT, but we want to make sure the peaks we see are in the correct positions, let us increase the size of the FFT by 5 times".
z = abs(fft(x,5*length(x)));
Linear convolution: In this case, you have to use a different size FFT. In linear convolution, a signal $x$ is convolved with a signal $y$. You might have heard that "Convolution in time is the same as multiplication in the frequency domain". This is true, but with the very important and often under-stated caveat that to do linear convolutions in the frequency domain, you have to make sure the FFT sizes are correct.

The correct size of the FFT lengths for a linear convolution, is $N_{\text{fft}} = L_x + L_y - 1$, where obviously those lengths correspond to the length of your signals $x$ and $y$. Why is this the case? Because in linear convolution, the result of convolving two signals is always equal to the sum of their lengths, minus 1. Do it on paper and you will see. Convolve a signal of length 3, and another of length 2. The result will be of length 4.

And so, as an example, one might say: "I want to linearly convolve two signals in the time domain, but I want to do this operation in the frequency-domain." They might do it like this:
Lx = length(x);
Ly = length(y);
Nfft = Lx + Ly - 1;
z = ifft ( fft(x,Nfft) .* fft(x,Nfft), Pfft);

Hope this helps.


Add odd/even harmonics to signal? How do I add odd or even harmonics to a floating point signal?

Do I have to use tanh or sin?

What I'm trying to do is achieve some very simple distortion effects, but I'm having a hard time finding exact references. What I'd like is something similar to what the Culture Vulture does by adding odd and even harmonics in its pentode and triode settings.  The float value is a single sample in a sample flow.

What your distortion box does is apply a non-linear transfer function to the signal: output = function(input) or y = f(x).  You're just applying the same function to every individual input sample to get the corresponding output sample.

When your input signal is a sine wave, a specific type of distortion is produced called harmonic distortion.  All of the new tones created by the distortion are perfect harmonics of the input signal:

If your transfer function has odd symmetry (can be rotated 180° about the origin), then it will produce only odd harmonics (1f, 3f, 5f, ...).  An example of a system with odd symmetry is a symmetrically-clipping amplifier.
If your transfer function has even symmetry (can be reflected across the Y axis), then the harmonics produced will only be even-order harmonics (0f, 2f, 4f, 6f, ...)  The fundamental 1f is an odd harmonic, and gets removed.  An example of a system with even symmetry is a full-wave rectifier.

So yes, if you want to add odd harmonics, put your signal through an odd-symmetric transfer function like y = tanh(x) or y = x^3.

If you want to add only even harmonics, put your signal through a transfer function that's even symmetric plus an identity function, to keep the original fundamental.  Something like y = x + x^4 or y = x + abs(x).  The x + keeps the fundamental that would otherwise be destroyed, while the x^4 is even-symmetric and produces only even harmonics (including DC, which you probably want to remove afterwards with a high-pass filter).

Even symmetry:

Transfer function with even symmetry:

Original signal in gray, with distorted signal in blue and spectrum of distorted signal showing only even harmonics and no fundamental:

Odd symmetry:

Transfer function with odd symmetry:

Original signal in gray, with distorted signal in blue and spectrum of distorted signal showing only odd harmonics, including fundamental:

Even symmetry + fundamental:

Transfer function with even symmetry plus identity function:

Original signal in gray, with distorted signal in blue and spectrum of distorted signal showing even harmonics plus fundamental:

This is what people are talking about when they say that a distortion box "adds odd harmonics", but it's not really accurate.  The problem is that harmonic distortion only exists for sine wave input.  Most people play instruments, not sine waves, so their input signal has multiple sine wave components.  In that case, you get intermodulation distortion, not harmonic distortion, and these rules about odd and even harmonics no longer apply.  For instance, applying a full-wave rectifier (even symmetry) to the following signals:

sine wave (fundamental odd harmonic only) → full-wave rectified sine (even harmonics only)
square wave (odd harmonics only) → DC (even 0th harmonic only)
sawtooth wave (odd and even harmonics) → triangle wave (odd harmonics only)
triangle wave (odd harmonics only) → 2× triangle wave (odd harmonics only)

So the output spectrum depends strongly on the input signal, not the distortion device, and whenever someone says "our amplifier/effect produces more-musical even-order harmonics", you should take it with a grain of salt.

(There is some truth to the claim that sounds with even harmonics are "more musical" than sounds with only odd harmonics, but these spectra aren't actually being produced here, as explained above, and this claim is only valid in the context of Western scales anyway.  Odd-harmonic sounds (square waves, clarinets, etc.) are more consonant on a Bohlen–Pierce musical scale based around the 3:1 ratio instead of the 2:1 octave.)

Another thing to remember is that digital non-linear processes can cause aliasing, which can be badly audible.  See Is there such a thing as band-limited non-linear distortion?


Relation between curve fitting and filtering Before I ask the question: I am from mechanical background and would appreciate a detailed answer.

Generally, for automobiles, I am trying to understand the following:

Given, say $1\;\mathrm{sec}$ of measurement data, which may contain a wide range of frequencies, we fit an $n^{\mathrm{th}}$ order polynomial curve through it. Doing this filters out certain frequencies which would have been present otherwise. It acts, as I can imagine,  like a low pass filter.

Question:

Let us say that we have a linear system with input $u$. Let $f$ be the input frequency. This linear system would produce an output, $y$, of same frequency, $f$. We fit an $n^{\mathrm{th}}$ order polynomial through these measurements. This curve-fitting acts as a low pass filter. What would be the characteristics of this low pass filter?

Context to my question and discussion:

How I see it is: The measurement data is $1\;\mathrm{s}$. Fitting a parabola through it can be approximated as a half of a sinusoidal curve. Similarly, a $3^{\mathrm{rd}}$ order curve would be a complete sinusoidal curve, and so on. So, for $2^{\mathrm{nd}}$ order polynomial, the cut-off frequency $f_c$ would be $0.5 \;\mathrm{Hz}$. Similarly, $1 \;\mathrm{Hz}$ for $3^{\mathrm{rd}}$ order polynomial and so on. Does this make any sense? Could somebody suggest me some book where this aspect is dealt with?

Generally speaking, I would not call a polynomial, in general,  a low-pass approximation to a signal. Per se, talking practically, standard definitions of the Fourier transform do not allow polynomials on $]-\infty,\infty[$ to have a useful spectrum. Generally speaking, polynomials (except constants) tend to infinity at fast rate, hence cannot be really limited in frequencies.

So it is difficult to say a polynomial filters out certain frequencies. However, locally, on an interval, they can pass nicely through data points, so they can smooth data, to some extend. Typical examples are sensor trends or drifts, when they are monotonous, and peaks or bumps.

You can find a lot of works on polynomial trend filtering, or polynomial peak filters. When the cost function for those fits is a squared error, you end up with a linear system in the polynomial coefficients, which turns into a weighted average of signal samples, that can be interpreted as coefficients of a linear filter. A famous example is the Savitsky-Golay filter, which was recently analyzed in frequency in 2011 by R. Shafer: What Is a Savitzky-Golay Filter?. This family of filters is also interesting as they provide smoothed derivatives.

But you can also fit a polynomial with other loss functions, like the $\ell_1$-norm, or with additional penalties. They are currently active topic in signal processing (I am working on some extension to Savitsky-Golay filters). Here are a few references, to start from:

On digital smoothing filters: A brief review of closed form solutions and two new filter approaches, 1986
Polynomial-Based Interpolation Filters---Part I: Filter Synthesis, 2007
$\ell_1$ Trend Filtering, 2009
Symmetrizing Smoothing Filters, 2013


What type of low-pass filter to use for gyroscope filtering in control application? I am using a gyroscope to generate a feedback signal that drives a motor to resist rotation, just as in a camera gimbal.

The gyroscope is sampling at 32 KHz with a bandwidth of 8 KHz, and my control loop is running at 1 KHz. Hence, I need to remove frequencies above 500 Hz. At the same time, I need to introduce as little delay as possible in order to be able to use a high control gain. I've attached the FFT of the sensor noise (source).

I am a little overwhelmed with the amount of choice available when it comes to filtering this signal. What filter properties should I be looking for? Are there any particularly suitable filters? This is implemented on an FPGA, so computational cost is not a constraint.

The lower the cut-off frequency of the anti-aliasing filter, the higher the delay thus degrading your phase margin. You need to be really careful when implementing an anti-aliasing filter in control loops applications.

First question, what is the amount of noise above 500 Hz? Based on the picture, you have noise at about 3 kHz, the main purpose of your anti-aliasing filter is to remove this spike, I would not worry too much about content between 500 Hz and 3 kHz.

Second question, what is the bandwidth of your control loop? 100 Hz? This is not the same as your sampling frequency. Typically it should be BW <= Fs/10. In your case BW <= 100 Hz.

Third question, what amplitude of aliasing can your application live with? If the bandwidth of your control loop is 100 Hz, aliasing occuring at f > 100 Hz will not really impact your performance, you only need to focus on aliasing that occuring at f <= 100 Hz. Even if there's aliasing occuring, if it less than say 1% of the signal of interest, you can probably neglect it.

Fourth question, can you run your control loop at 32 kHz? You would not need an anti-aliasing filter then. (I'm assuming that there's already an analog anti-aliasing filter in the gyroscope).

Based on the information we have, here's what I would do :

If you cannot sample faster than 1 khz, the biggest problem is that the 3-kHz spike will alias to DC (or near DC) . So you need to filter it. I would design a notch filter that removes your 3 kHz spike. This notch filter will add a smaller delay compared to a low-pass filter with a 500 Hz cut-off.

However, if you can sample faster, say 4 kHz, then the 3 khz spike will alias at 1 kHz which should be outside the bandwidth of your controller. In that case, I would only design a simple low-pass filter with a high cut-off frequency.

Edit : You mention you have an FPGA. In that case keep your sampling at 32 kHz and design a controller with a Bandwidth of like 100 Hz. In that case you don't need to worry about aliasing and noise since the noise will be outside the bandwidth of the control loop.


Determining type of filter given its equation Given a general filter equation, how can one tell the type of filter that the same equation represents? Meaning how can I know if the filter is Low/High/Band Pass etc...?

For example, the following equation:

$ y_{n} = x_{n} + x_{n-2} $

Represents a "Band Stop" filter, but why?

The $\mathcal Z$-transform of the given difference equation is :

$$H(z) = 1 + z^{-2}$$

Put $z= \exp(j\omega)$ and remember that in discrete time systems the low frequencies are at $2n\pi$ $(n=0,1,2, \ldots)$ and the high frequencies occur at $(2n+1)\pi$; this is just a consequence of the periodic behavior of the discrete time complex exponentials.

So, $H(z) = 1 + \exp(-2j\omega)$ at $z = \exp(j\omega)$.

When $\omega=0; H(z) = 2$ and $w = \pi$ gives $H(z)= 2$.

Thus, both at high and low frequencies the the system function provides same
gain and hence the filter with the given $H(z)$ is a BAND REJECT/ NOTCH FILTER with $H(z) = 0$ at $\omega = \pi/2$.


Basic Confusion About the DFT and Convolution I am learning DSP (with Digital Images) and I have some elementary confusion about the convolution between two discrete periodic signals. Specifically, I have learnt that when filtering an image, we apply convolution between the image (as discrete 2D signal) and a kernel filter. To visualize what happens in the frequency domain, we calculate the DFT of the image and the kernel and since we are convolving in the space domain, in the frequency domain this corresponds to a point per point multiplication of the two spectrums.

Now, the DFT assumes that the signal to transform is discrete AND periodic, resulting in a spectrum which is discrete and periodic. So here we are assuming that both the signal and the filter are discrete and periodic signals which repeat (usually) after the same M samples. My confusion is, how can convolution work for two discrete periodic signals? Since they repeat infinitely, applying the convolution usual definition as a flipped-shifting multiplication and sum, we would end up with an infinite sum in some single points if the signal to be filtered and the filter overlap. (Because they would overlap in the current repetition, in the next, and the next ad infinitum)

As an example in 1D, let $x[n]$ be our discrete function, which is periodic with a period of $N$ samples in the spatial domain. And let $h[n]$ be the filter impulse response in the spatial domain, also periodic of $N$ samples. We want to filter $x[n]$ with $h[n]$ and get the filtered signal $y[n]$. To do this we convolve these two periodic discrete signals:

$$y[n] = \sum_{k = -\infty}^\infty x[k]h[n - k]$$

However this convolution for a certain $n$, e.g. $n=0$, can lead to an $y[0]$ which is infinite, as both $x[n]$ and $h[n]$ are periodic.

When we convolve between the image and the kernel, we are only considering one repetition of $h[n]$ and one repetition of $x[n]$. Like this:

$$y[n] = \sum_{k = 0}^{N-1} x[k]h[n - k]$$

Why can we do this? Will this still give us a periodic discrete $y[n]$ whose representation in the frequency domain is the point-per-point multiplication of the original signal and filter?

As I understand your question, it boils down to: why is it ok, to take the convolution of a periodic signal only on an interval rather than over infinity? The DFT or discrete nature of the signal does not really play a role in your question. It merely provides the basis of assumung your two signals $f$ and $g$ as periodic, on the interval $[0,T]$ such that $$f(t+T)=f(t-T)=f(t)$$ and obviously $g(t+T)=g(t-T)=g(t)$.

Now assume that you define the convolution over the whole $(-\infty,\infty)$ like so:
$$y(t):=(f\ast g)(t)=\int_{-\infty}^\infty d\tau\, f(t-\tau)g(\tau)$$

From the periodicity above follows that
$$\int_{-\infty}^\infty d\tau\, f(t-\tau)g(\tau) = \int_{-\infty}^\infty d\tau\, f(t+T-\tau)g(\tau)$$
which you can intepret intuitevely to mean that the convolution should also be periodic. I write "intuitively" because the integral is divergent as we will see below.

Furthermore you can split the integral into a sum of integrals where $\tau$ goes over intervals of length $T$, such as $[n\cdot T,(n+1)\cdot T]$ with $n$ a signed integer number in $(-\infty,\infty)$. Doing that we have (after some simplifications)

$$\begin{eqnarray}
y(t)&=&(f\ast g)(t)=\int_{-\infty}^\infty d\tau\, f(t-\tau)g(\tau) \\
&=&\sum_{n=-\infty}^\infty \int_0^Td\tau\,f(t-\tau)g(\tau)
\end{eqnarray}$$
which would be divergent just because of the sum. Thus, the convolution of periodic functions produces a divergent integral. Sticking to this definition of the convolution would make it impossible to define the convolution like that. At least one of the functions must go to zero sufficiently fast, to produce a finite integral. Your periodic functios, however, will violate this decay condition. Thus it makes sense to window the functions with a rect function that captures all the information in one period.

So you see why it would make sense to define the cyclic convolution as just the integral

$$y(t)=\int_0^Td\tau\,f(t-\tau)g(\tau)$$

and use the periodic extension, since we already provided some intuition as to why setting $y(t+T)=y(t)$ makes sense. Thus the cyclic convolution carries all the information you would have from the regular definition but avoids the divergence due to the infinite sum. I hope this helps.


Is there any good book regarding Image Processing Algorithms for Feature Extraction and Matching which gives examples in MATLAB? I have been looking to implement feature extraction and matching in MATLAB. For this, I've read many algorithms, but i'm not able to understand how to implement these in MATLAB. Is there any Book related to this topic which is based on MATLAB? i.e. gives exapmles and sample codes in MATLAB?
I have read Gonzolez and there is no other book that deals with Image Processing Algorithms with help of MATLAB.

you may refer Feature Extraction & Image Processing
 By Mark Nixon, Alberto S Aguado.


What's the Q factor of a digital filter's pole? For the analog S-plane, the Q factor depends on the angle of the complex pole $p$ from the horizontal axis, $\theta = \arg p$:

damping factor $\zeta = \frac{1}{2Q}$, and $\zeta = \cos(\theta)$, and $\theta = \arg(p)$ so $Q = \frac{1}{2 \cos(\theta)}= \frac{1}{2 \cos(\arg(p))}$  Right?

But in the Z plane, the curves of constant Q are more complicated:

So how do you get the Q of a complex pole?  To transform an analog pole to the Z plane, you use $z = e^{sT}$?  But when I try to apply this I must be doing something wrong.

Oh, maybe I got it.  Maxima was making assumptions about the pole location being real?

$$Q=-\frac{\sqrt{{\log(|p|)}^2+{\arg(p)}^2}}{2\log(|p|) }$$

it depends on how you map the analog filter to the digital filter and how the s-plane poles get mapped to the z-plane poles (ya know, the "bilinear transform" vs. "impulse invariant" vs. whatever else is out there).

probably, if it were up to me, for the purpose of defining $Q$, i would map the poles as you would map $s$ to $z$ with

$$ z = e^{sT} $$

where $T \triangleq \frac{1}{f_\text{s}}$ is the sampling period and $f_\text{s}$ is the sample rate.  that would be assuming impulse invariant.  the damping rate of the digital filter would be the same as the analog filter because both have the same exponential impulse responses (it's just that the digital filter's impulse response is the sampled analog impulse response).

for bilinear transform the $s$ to $z$ mapping is

$$ z = \frac{1 + sT/2}{1 - sT/2} $$

if you're using $Q$ to define the degree of resonance (like the dB boost you get around the resonant frequency, whether the filter is LPF, BPF or HPF), then this mapping will preserve the height of that bump for a given $Q$.

i dunno if another mapping would be useful for this definition of digital Q.  in the Audio EQ Cookbook, i worried a little about mapping BPF bandwidth from s-plane to z-plane and tossed in a fudge factor because bilinear transform sorta scrunches the bandwidth as the resonant frequency gets closer to Nyquist, $f_\text{s}/2$.  if preservation of bandwidth, as you map from s-plane to z-plane, is important to you, then i would take a look at that.

so, given the s-plane pole location (which is a function of $Q$), map it to the z-plane and you'll get an expression for $Q$ based on the real and imaginary parts (or magnitude and angle) of the z-plane poles.


Kalman filter with multiple sensors I have been reading more about EKF and I am a bit confused on how you handle predictions with multiple sensors.
Ex, I have IMU, GPS, Odom and Stereo Camera.
Each can be used to predict location, how do you combine them?
My understanding of EKF fits in the basic framework of predict, act, measure, fuse prediction with sensor.
But it's not clear how multiple sensors that measure the same thing is used.

Summary
You can have only one (or no) primary sensor which is used in the prediction process. All other sensors are secondary and only used for corrections. Each correction takes  the current state and a measurement from one of the secondary sensors and outputs the corrected state. You can run corrections as often as you want, typically each time one of the secondary sensors produces a new measurement. There is no need to run corrections at a fixed interval.
The new state is predicted from the model matrices F and B and from the old state x and primary sensor data (or control input) u.
$$x_{k+1} = Fx_k + Bu_k$$
The corrected state is calculated from the kalman gain matrix K, observation matrix H and secondary sensor data z. Each secondary sensor has its own observation matrix.
$$x_{k,corrected} = x_k + K_k(z_k-H_k x_k) $$
Detailed version
You have to designate one primary sensor each for orientation and position. Alternatively, you can go without primary sensors if system dynamics allow. The primary sensors should be extremely reliable, accurate in the short term, and support a high update rate. Long term drift is not a concern here.
In a car navigation system, the primary sensors are usually steering angle and odometer, the secondary sensor is GPS. In an attitude and heading indicator for aircraft, the primary sensor is a three-axis gyro and secondary sensors are accelerometer and magnetometer.
There can be an arbitrary number of secondary sensors. Their purpose is to provide long term accuracy, which the primary sensors lack. The secondary sensors have other shortcomings, like unreliability (no GPS in tunnels), poor short term accuracy (apparent gravity vector disturbed by g-forces from maneuvering), or incomplete information (gravity vector provides pitch and bank, but not heading). If you have a sensor that can take both roles - a sensor that is reliable, provides enough information to uniquely deduce orientation or position, and is accurate both in the short and long term, then you don't need a Kalman filter in the first place.
The prediction stage of the Kalman filter runs at a fixed and high rate and continuously updates the current estimates for position and orientation based on the old estimates and measurements from the primary sensors. If you operate the Kalman filter without primary sensors, it turns into a simple model-based extrapolator. Needless to say, extrapolating from the last known GPS position using steering and odometer data is much more accurate than a simple constant-speed extrapolation.
The correction stage of the Kalman filter becomes active whenever a measurement from a secondary sensor arrives. Each correction handles the measurement from only one sensor and updates the current estimates for position and orientation. There can be an arbitrary number of corrections with arbitrary or even irregular update rates (still no GPS in tunnels, therefore no corrections for several minutes).


Discrete Time Fourier Transform (DTFT) for an unstable system (Ideal Low Pass Filter) The Dirchlet conditions state that if the signal is absolutely summable then it the DTFT of the signal definitely exists. This is a sufficient condition but not necessary condition.

There are systems like Ideal Low Pass Filter, which are not absolutely summable. The condition of absolutely summability is a necessary and sufficient condition for stability. Therefore, I can find a bounded input for which the system produces unbounded input.

For such systems, which are not absolutely summable, but are absolutely square summable, we define the DTFT in the same manner and argue that even though non-zero error exists, the energy of the error is 0, which essentially gives rise to Gibbs Phenomenon.

But the condition of $y[n] = x[n] * h[n]$ is independent of stability, and simply depends on Linear and shift invariance property. From this definition, I can find a bounded input which will produce unbounded output. But the relation, $Y(w) = X(w)H(w)$ is also independent of stability. Even if I take into account Gibbs Phenomenon, I am unable to deduce the existence of a bounded input which will produce and unbounded output, and this essentially suggests the fact that bounded input will give bounded output.

Where am I going wrong?

It is generally not true that the relation $Y(\omega)=X(\omega)H(\omega)$ is independent of the system's stability. For systems with a rational transfer function (i.e., systems that can be described by differential or difference equations), the frequency response, which is the Fourier transform of the impulse response, only exists for stable systems. For marginally stable systems with poles on the imaginary axis (for continuous-time systems), or on the unit circle (for discrete-time systems), the Fourier transform exists if we allow Dirac delta impulses in the expression for the Fourier transform. So for marginally stable systems you can see that the output can become unbounded, even for bounded input signals, due to the Dirac delta impulses in their frequency response. For unstable causal systems with transfer function poles in the right half plane, the frequency response doesn't exist.

For systems for which there exists no transfer function, but which do have a frequency response, such as ideal frequency-selective filters (low pass, high pass, etc.), the relation $Y(\omega)=X(\omega)H(\omega)$ holds, regardless of stability. Ideal frequency-selective filters are unstable (and non-causal), but their frequency response $H(\omega)$ is finite. This means that the frequency response of the bounded input signal which generates an unbounded output (if it exists) must be unbounded for at least one value of $\omega$, and the inverse Fourier transform of $X(\omega)H(\omega)$ must result in an unbounded time-domain signal. Also note that there is no guarantee that the Fourier transform $X(\omega)$ of the bounded input signal that generates an unbounded output signal even exists. And if it exists, it's generally not immediately clear from the expression $X(\omega)H(\omega)$ that the corresponding inverse transform is unbounded.

I think it is a misunderstanding that the relation $Y(\omega)=X(\omega)H(\omega)$ somehow implies that a bounded input must result in a bounded output. It's just not so easy to see whether the output is unbounded because one would have to compute

$$\frac{1}{2\pi}\int_{-\pi}^{\pi}X(\omega)H(\omega)e^{jn\omega}d\omega$$

and check for each value of $n$ if the expression becomes unbounded. Direct time domain analysis is more straightforward in this case.


About understanding a Matched Filter I have some basic doubt in understanding a matched filter (MF). Please note that I have already seen very informative post in this context Understanding the Matched Filter

In what follows, I list three information about matched filter that I found in standard textbook and wikipedia:

(1) Wiki page says that MF is obtained by correlating a known signal, or template, with an unknown signal to detect the presence of the template in the unknown signal.
In digital communication, I assume that the known template is the pulse shape p(t) used by the transmitter [e.g., 'sinc-like' impulse response of root-raised cosine (RRC) spectrum] and the unknown signal is the noisy received signal r(t).
Assuming that the receiver knows p(t) used by the transmitter, the receiver convolves r(t) with p*(-t). The output of convolution is sampled at decision points.

(2) A standard textbook on communication says that the receiver does a convolution of the received waveform with the time-reverse conjugate of channel response h(t). That means, MF operation consists of convolution of r(t) with h*(-t), followed by a sampler.
(3) The textbook also says that MF can be seen as an equalizer.

I understand that MF does a convolution followed by a sampling. That is it. Simple. But I am seeking detailed clarity in my understanding of the philosophy of this operation.

Questions:

(A). What is p(t) in the MF operation? Suppose I have NRZ data, then do I use the 'rectangular pulse of NRZ' or the 'sinc-like impulse response of RRC' as p(t) in the convolution? Wiki says, rectangular NRZ pulse is used. Why so?

(B). Refering to the two different notions in (1) and (2), seems my understanding of the relation between p(t) and h(t) is unclear. Can someone make it clear with simple explanation?

(C). How do you explain MF operation as an equalization? Equalizer removes ISI, MF removes the effect of noise...

(D). What is partial MF?

I really appreciate your help in improving my understanding of this beautiful concept of matched filtering.

I just want to comment on the confusion between $p(t)$ and $h(t)$. What happens is that we transmit a signal $s(t)$ where

$$s(t)=\sum_ks_kg(t-kT_s)$$

where $s_k$ are the data symbols, $g(t)$ is the transmit pulse shape, and $T_s$ is the symbol time.

Transmitting this over a wireless channel with a channel impulse response $h(t)$, the channel acts acts as a filter, and thus the channel output is the convolution between the transmitted signal and the channel, i.e.,

$$y(t) = s(t)\star h(t) = \sum_ks_k\underbrace{\left[g(t-kT_s)\star h(t)\right]}_{p(t-kT_s)}$$

So, the convolution of the transmit pulse shape and the channel is considered in the matched filter design.

By the way, the purpose of the matched filter is not to "remove the noise", but rather to maximize the signal-to-noise ratio (SNR) by maximizing the received signal power at the output of the matched filter.


Isosurfaces from three dimensional column data: methods I have just been asked the following question, and   I somehow felt short of smart answers. You are given a series of $N$ triplets of values ($P_1$, $P_2$, $P_3$), pertaining to physical measurements. $N$ is in the order of $100$ to several thousands. My colleague would  like to draw isosurfaces from this dataset.

I have not more (yet) information about the quality of space filling, and the uncertainty around the data.
So I suggested to interpolate the given point cloud to a regular grid and use standard isosurface extraction. But there might be other, more direct methods, dealing with point cloud issues.

Some related papers:

Direct Isosurface Extraction from Scattered Volume Data
Isosurface Generation for Large-Scale Scattered Data Visualization

Would you suggest sounder approaches (references, software implementing them are a plus)

You are given a series of $N$ triplets of values ($P_1$, $P_2$, $P_3$), pertaining to physical measurements. $N$ is in the order of $100$ to several thousands. My colleague would  like to draw isosurfaces from this dataset.

This is a little bit broad but perhaps still with some chances of accurate answers.

To an extent, it is important to know where the data comes from because it affects how the iso-surfaces are reconstructed. A Digital Elevation Model, for example, is essentially a 2D construction but volumetric data are inherently 3D. In the latter case, a trully 3D surface has to be fitted through the points of the contour which is much more challenging.

In both cases, the marching cubes algorithm is the standard way to extract the iso-surface but other methods involving a shrinking elastic surface (basically, 3D generalisations of active contours) have also been used.

In terms of practically applying it:

For relatively small point-clouds (a few hundreds of points), Blender with Point Cloud Skinner, can produce quick results. I have used this to create a surface through Magnetoencephalography sensors (whose locations are known) which ends up looking like a helmet. If your point cloud describes a convex shape, this will have no problems recovering it. I think that it works similarly to QHull (i.e, extracts the convex hull) but it has a threshold parameter that to an extent allows it to account for holes.
Paraview. Incredibly useful and with Python scripting capabilities. The software handles loading the data and several standard operations and allows extracting the results. The good thing about Paraview is that it includes a "raw import" functionality which you can use to load data in whatever format you may have it in. Isosurfaces are produced with "Contour". This is straightforward for volumetric data. If you already have the points, there is a "Table to structural grid" function that will bring the data to the expected form prior to applying "Contours".

I would consider this the faster option over the Visualisation Toolkit but if you need more control over the result, you might find the VTK more useful.

Meshlab. Meshlab is specifically built for handling huge point clouds including fitting surfaces. Meshlab is interesting because there are various scripts you can use around it, depending on what you are trying to achieve.
The Point Cloud Library and its Python bindings. This is similar to Meshlab in that it can handle huge volumes of data and it includes a surface fitting capability. It is however more challenging in setting up and operating.

Hope this helps.


Algorithms to re-spatialize a stereo recording audio signal? Let's say we have an audio signal stereo recording, made for example with a XY microphone positioning:

Which algorithms can be used to "re-spatialize" this recording, i.e. try to virtually "move the microphones", and recreate a new stereo signal, for example with an AB mic positioning?

To do this, we would need to be able to:

know the effect on the signal of rotating a microphone
know the effect on the signal of moving a microphone a few centimenters from each other

I can imagine that this involves delaying/shifting x_L[n], x_R[n] (phase), but more generally should we use convolution / STFT techniques?

Which algorithms can be used to "re-spatialize" this recording, i.e. try to virtually "move the microphones", and recreate a new stereo signal, for example with an AB mic positioning?

In general, this is a beamforming "problem", but it cannot be done in exactly the way it is described here.
With a setup like this, you can adjust the relative differences between signals recorded on each microphone, but you would not be able to tell if rotating the microphones (or putting them apart) would result in one of the microphones:

Moving closer/further away from a source.
Moving closer/further away from a local pressure maximum/minimum in the room.

The way to overcome 1 (and partially 2) is to use an array of $M$ microphones and beamforming techniques. You would first localise the sources and the signals they contribute to the array and then place a set of $N<M$ microphones within the array, re-project the signals and see what this sounds like.
This deals with 1 and it deals with 2 as long as, the microphones (or sources) do not move too far away from their locations.
A microphone does not only pick up the sound that is produced by a source in its vicinity but it also picks up the interference of the space the sound propagates in on that sound.
As the sound waves propagate, they get reflected, refracted, absorbed over surfaces and objects in a room. This is what gives the room a "signature". You can tell for example a recording done in an empty warehouse because it has more reverberation. You can tell if the source was near or far from the microphone because in the first case the reverberation is softer and it is trailing the main sound but in the second case, the reverberation is thicker and the main sound has receded deeper into its reverberations.
If you tried to use beamforming to record sound in a tunnel and then asked, "well what would it sound like if we placed the microphones further back in the tunnel?" (many more reflections) Then, it is impossible to work out the dimensions of the tunnel and the effect these would have on the sound the mics pick up (without assumptions or some other way of filling your knowledge gap).
You would "rotate the phase" (for the new configuration of the mics) but you would not know the "phase correction" required at a given point in space because of the shape of this space.
For more elaborate examples of this type of beamforming, please see here or here and more generally look for work around audio beamforming.
Hope this helps.
EDIT:
On the points made more specific in the comments:

Are there formulas about how to shift signals according to a distance change?

Yes, that would be the typical way of working out "phase of arrival" taking into account the speed of sound in air and frequency of a component (from which you workout the wavelength and through that the frequency difference because of distance). These equations are present in the beamforming examples.


How to determine stopband of discrete Gaussian, stdev sigma, support N I would like to analyse a gaussian smoothing kernel with a set standard deviation and support (let's say, in MATLAB, fspecial('gaussian', [5 1], 1.3) so sigma is 1.3 and support is 5) in the DTFT domain to determine the stopband.

My understanding is that the stopband is when the magnitude of the impulse response is below -3 dB. However, I don't have a clue how to determine this given a sigma and support value.

I'm particularly interested to see if the DTFT quantifies the artifact due to the compact support as I would like to understand this a little better.

I found answer at https://courseware.ee.calpoly.edu/~fdepiero/curr/dsp/dsp11.htm . Converted to the present problem, in MATLAB:

kernel2 = fspecial('gaussian', [1 5], 1.3)
spike = [1 0 0 0 0];
[h, w] = freqz(kernel2, spike, 1024);
m = abs(h)
m_db = 10*log10(m/max(m));
plot(w,m_db)

Because the compact support Gaussian had zero crossings, I used the first zero to define the stop band.


What are the effects on filtering upsampled data without interpolation? I have an application that displays several signals from different sensors at various sample rates.  In order to display the data, I "stretch" the signal by repeating samples to match the highest sample rate.  For example: if I display a 500 Hz signal and a 100 Hz signal at the same time, every sample in the 100 Hz signal will be repeated 5 times so it appears to have a 500 Hz sample rate.

When I apply filters to the signals (all second-order IIRs: LPF, HPF, notch) I use filters designed for the highest sample rate.  In my example here, I'm using 500 Hz filters on an upsampled 100 Hz signal.

How "bad" is this?  Is my only option to design an interpolation filter during upsampling, and if so, can anyone lead me to good references on how to design such a filter?

How "bad" is this?

It's not too good.  What you are doing is equivalent to upsampling by a factor of five (i.e. inserting four zeros in between every sample), and then filtering it with a filter that has the impulse response [1, 1, 1, 1, 1].  The picture below shows the frequency response of this filter.

The signal aliases introduced by the upsampling are centered at 0.4 and 0.8, so the good news is that the filter has a null at the center of the aliases.  The bad news is that unless the signal has a really narrow bandwidth, the signal aliases will be much wider than the nulls, and a lot of alias energy will distort your signal.

  Is my only option to design an interpolation filter during upsampling, and if so, can anyone lead me to good references on how to design such a filter?

If you want to understand how to design a low-pass filter then I would suggest learning about windowed sinc functions.  If you just want to get a filter that works then I would recommend installing Octave and using the fir1 function to design the filter.


Why would one use a Hann or Bartlett window? Suppose we're designing a low-pass FIR filter, and I want to use one of these three windows: Bartlett, Hann or Hamming. From Oppenheim & Schafer's Discrete-Time Signal Processing, 2nd Ed, p. 471:}

All three of them provide the same transition band width: $$\Delta\omega=\frac{8\pi}{N}$$ where $N$ is the order of the filter and is assumed large enough.

However, the overshoot (let's call it $\delta$) is different for each window, an the following inequality holds:

$$\delta_{Hamming}<\delta_{Hann}<\delta_{Bartlett}$$

So if we use a Hamming window, we get the smallest overshoot and a transition band with width $\Delta\omega$. If we use one of the other two windows, the width of the transition band is the same, but the overshoot increases.

This leads me to think that there is no case in which one would use a Hann or a Bartlett window, as the Hamming one is better than them: it improves one aspect ($\delta$), stays the same in another ($\Delta\omega$).

The question is: why would someone choose a Hann or a Bartlett window if a Hamming one can always be used?

In reviewing Fred Harris Figures of Merit for various windows (Table 1 in this link) the Hamming is compared to the Hanning (Hann) at various values of $\alpha$ and from that it is clear that the Hann would provide greater stopband rejection (The classic Hann is with $\alpha =2$ and from the table the side-lobe fall-off is -18 dB per octave). I provided the link as you can see many more considerations involved when choosing a window for various applications.
The result of this is apparent when comparing the kernels for a 51 sample Hann and Hamming window using Matlab/Octave. Note the higher first sidelobe level with Hann but significantly greater rejection overall:

Personally, I would not use either window for filter design. If any window, I would use the Kaiser window, or preferably firls. See FIR Filter Design: Window vs Parks McClellan and Least Squares for the related discussion.
I convolved a 26 sample Hann with a 26 Hamming to come up with an alternate 51 sample "Hann-Hamming" with the following result:

UPDATE: This Hann-Hamming does not (generally) out-perform a Kaiser window of similar main-lobe width:

I then tried what I call a "SuperKaiser" where I convolved two shorter length Kaiser windows to come up with an alternate 51 tap window with the following result. This was done by convolving Kaiser(26,5.5) with Kaiser(26,5.5) such that SuperKaiser(51,5.5)= conv(kaiser(26,5.5),kaiser(26,5.5). At first glance it appears to generally outperform the kaiser(51,12), matching the main lobe width and providing superior stopband rejection over most of the stopband. An integration of total stopband noise under the assumption of AWGN is of interest to see if this new window is superior under that condition (does the relative area under the first two sidelobes where SuperKaiser is inferior completely offset all the remaining stopband improvement?). If I have time I will add that assessment. Interesting! As @A Concerned Citizen astutely pointed out, as we keep convolving we should approach a Gaussian window.


FFT analysis of binned data of different length I am conducting experiments to collect wind speed data from wind anemometers placed on a moving platform. Closely, there is a fixed wind mast holding a wind vane. Prior to analysing the wind data from the moving platform, I binned the data with 1-minute average wind direction. Each bin contains a number of continuous data sets (Set 1, Set 2, Set 3 etc.) which may not necessarily have the same number of data points. For example, Set 1 has 20,000 data points, Set 2 has 6,000 data points, etc.

I am performing an FFT on each of the sets to detect any periodicities in the signal. I am plotting the real amplitude of the sinusoids present in the signal on the y-axis vs. the frequency of the sinusoid. The FFT of each data set is yielding different amplitude results, on the y-axis and I'm not sure whether this is a result of the different length of the arrays. I am interested in an average value of the y-axis of the spectrum. However, I don't think that it is suitable to just calculate the average due to the different array size. Can someone please suggest a method? I am not used to analysis of signals in this way so I'm a bit confused. Thanks a lot.

Let's say you have several data sets $x_i$ and $N_i$ is the length of the $i$th data set.

First you're right that the amplitude of the FFT output scales with its length. To normalize all data sets they have to be divided by $N_i$. (Note that depending on the implementation of the FFT there might be a scaling involved already, check the Doc or see this question)

To combine all data sets you should append $(\max_i N_i) - N_i$ zeros to every data set before taking the FFT. The frequency axes are then aligned correctly and you can calculate the average of all FFT outputs.

Btw, I would plot the absolute value of the FFT output, not the real part. And sometimes I find it useful to "smooth" the spectrum by averaging every $M$ neighbouring FFT bins. Last, I have assumed that the data of every set has been acquired at equidistant time steps, otherwise the FFT isn't applicable at all.


How to find the normalized coefficients? I use the below matlab function to design the unquantised coefficient set:

  [b,a]=ellip (4,2,60,2*[0.2,0.3]);

and then i use the matlab finction "tf2sos" to convert the direct form coefficients to a set of a factorized second-order sections:

  [sos,g]=tf2sos(b,a);

My result is:

The normalized version of the difference equation is:

How i can find the b and the c1,c2,...

Is the actual result in this form:

or factorized means normalized and the result is in this form:

The gain plays any special role?

I'm not sure how you define your normalized difference equation, but in any case you can equate the $b_{ij}$ coefficents in the sos matrix with your $c$ coefficients, and the gain that tf2sos returns would then be your normalization constant $b$.

I think you already know, but just to be clear: every row of the sos matrix contains the coefficients for one second-order section, starting with the 3 numerator coefficients, and then the denominator coefficients (the first of which is always 1). You can incorporate the gain into the first section.


Why do you get instantaneous phase by calculating correlation with a carrier signal? I am somewhat new to DSP. I'm trying to implement BPSK by following this Python QAM tutorial.

My code seems to work, but I really just ported the tutorial's Python code to C++, and do not understand it. This is the part that confuses me the most:

# Demodulate in "real" and "imaginary" parts. The "real" part
# is the correlation with carrier. The "imaginary" part is the
# correlation with carrier offset 90 degrees. Having both values
# allows us to find QAM signal's phase.

for sample in samples:
    t += 1
    angle = CARRIER * t * 2 * math.pi / SAMPLE_RATE
    real.append(sample * math.cos(angle))
    imag.append(sample * -math.sin(angle))

for t in range(0, len(real)):
    # This converts (real, imag) to an angle
    angle = math.atan2(imag[t], real[t])
    phase.append(angle)

I understand that this code is correlating the input signal with two identical carrier waves which are offset by 90º. I do not understand why that correlation has anything to do with phase.

I'd love an answer to any of the following questions:

Why does this code work?
What is the connection between correlation and phase?
Is this a common, reliable method of calculating phase?
Why does the FFT need so many samples to calculate phase, but this code just needs 1 sample?

Feel free to refer me to a signal theory or DSP book that answers these questions!

This is the analytic signal. See wikipedia for better details.

Basically you are taking the carrier correlation, $f(x)$, and its Hilbert transform, $(h \ast f)(x)$ and combining them into a single complex valued signal, the analytic signal

$$f_A(x) = f(x) - i (h \ast f)(x)$$

which can be written as the complex exponential

$$f_A(x) = A(x) e^{i \phi(x)}$$

Why this is useful is that it performs a split of identity into seperate amplitude $A$ and phase $\phi$ components. Actually the Hilbert transform requires the entire signal to compute, but because you know the function of your signal already you only need the one sample.To localise the phase, you would use a bandpass filter.


How to Mesure the smoothness of a signal so how can i determine if a signal is smooth or not ? And if its possible to get something indicating the level of the smoothness of my signal.

I looked at the r-squared and hurst exponent but i wanted to know if they where others methods.

Hurst exponent is a good option.

For a better answer you must better explain your objective. For example if you are comparing signals and if they have the save length, amplitude, sampling frequency and/or others.

If you are looking for very simple approaches then, for example, you can:

calculate the discrete derivative
calculate absolute value
sum all values

Another option could be:

apply a low pass filter
calculate difference between original signal and filtered signal
a smaller difference means a smoother signal.


Why does the solution for the voltage of this RLC circuit result in a maximum capacitor value? I have a circuit with an inductor (L), capacitor (C), and two resistors (R1=R2):

This represents an equivalent electrical circuit to a plectrum or hammer (cap=spring, inductor=mass) striking a string divided into two segments (R1 & R2).

The system is energized by a current impulse at time 0. This DSP synthesis technique is explained more here: https://www.dsprelated.com/freebooks/pasp/String_Excitation.html

The problem I am having is that the solution for the voltage leads to a condition where there is a maximum capacitor value or it is unsolvable, and I am not sure why or what this represents.

The current across the parallel component (the "exciter") is described as follows:

$I(s) = I_C(s)+I_L(s)$

$I(s) = sCV_C(s) + \frac{1}{sL}V_L(s) + \frac{i_0}{s}$

$I(s) = \frac{CLs^2+1}{Ls} V_{exciter}(s) + \frac{i_0}{s}$

The current across each string is represented by:

$I(s) = V_R(s)/R$

Furthermore, if $R1=R2$, $V_R(s)$ can be approximated given that with the exception of the momentary impulse:

$V_{exciter}(s) + V_{R1}(s) + V_{R2}(s) = 0$

$V_{exciter}(s) + 2V_{R}(s) = 0$

$V_R(s) = \frac{-V_{exciter}(s)}{2}$

Therefore we can put these equations together as:

$\frac{-V_{exciter}(s)}{2R} = \frac{CLs^2+1}{Ls} V_{exciter}(s) + \frac{i_0}{s}$

$V_{exciter}(s) = \frac{-2i_0LR}{2CLRs^2+ Ls + 2R}$

To then convert this to the time domain, I need it in terms of: $\frac{1}{(s+a)(s+b)} = \frac{e^{-at}-e^{-bt}}{b-a}$

This is solved as:

$V_{exciter}(s) = \frac{-2i_0LR}{2CLR(s^2+s/(2CR)+1/(CL))}$

$V_{exciter}(s) = \frac{-i_0}{C}\frac{1}{(s+a)(s+b)}$

$a=\frac{-1}{4LR}\left(-\frac{L}{C}-\sqrt{\frac{L^{2}}{C^{2}}-\frac{16CLR^{2}}{C^2}}\right)$

$b=\frac{-1}{4LR}\left(-\frac{L}{C}+\sqrt{\frac{L^{2}}{C^{2}}-\frac{16LR^{2}}{C}}\right)$

This results in a final solution of:

$v_{exciter}(t) =\frac{-i_0}{C}\frac{e^{-at}-e^{-bt}}{b-a}$

The problem now is the square root term for $a$ and $b$. This creates a condition where in order to be solvable:

$L (L - 16 C R^2) >= 0$

So at a given $L$ and $R$ value, $C$ must be $< \frac{m}{16R^2}$ for this equation to be solvable or you otherwise end up with a square root of a negative which I am not sure how to deal with and get real results from.

Why does this happen and what does it mean? Is there some way to solve this equation differently so it will still work at any $C$ value?

It was suggested to me that above this capacitor threshold, the circuit may be going from exponential decay to sinusoidal.  But even still, should I not be able to work out an equation for the voltage that will still give me a value?

Why the fear of complex numbers? Just do a distinction of cases and generate solutions for the real and the complex case. $b-a$ in the denominator will always be real, as the imaginary parts will cancel out each other, so complex numbers will only be in the exponents, wich will lead to some nice sinusoidal functions, which is expected, as this is an oscillating system. The real case just represents the overdamped case with no oscillation. If you don't know what I'm talking about have a look at this wikipedia article:
https://en.wikipedia.org/wiki/Damping_ratio


Measuring time distance between bass tones and send it as a signal I want to measure in my file the time distances between two or more bass tones.

Explanation:

I took a YouTube video and converted it to a .wav file.

As you can already see in the title, the bass was manipulated, before it was uploaded. (Probably with audacity or something else). I cut out about 16 seconds, from the beginning of the video and read the audio file in matlab. Now in theory, (i believe so) bass (lower frequency) generally have a higher amplitude compared to the other tones. A higher amplitude means, that the tone is louder. So basically, if i do a fft on it, i get the frequency spectrogram of my signal, which is represented in the time domain. Then i square the frequency spectrogram with 2. I dont know, how is it called, since english is not my mother language, but I believe this is the absolute value in the mathematics or also the power from the frequency spectrogram. As i mentioned, the amplitudes of bass tones are very high, therefore i should be able to see, when a bass tones occurs in the spectrogram.

I want to measure the time distance between two or three bass tones. If i get this, then i want so send this values in the form of signal through an Arduino Uno, which controls a vibration motor. If the bass comes, for example, every two seconds, then the vibrations motor will/should vibrate every two seconds too.

Question:

I want to measure the time distance between two bass notes and send it as a signal to my motor. Depending on the music type of my file. What are the necessary steps, to do that? How do i have to write my code line? I have already tried some possiblities, but i don't know further and i want to realise it with matlab, since i am not experienced in other languages. I can alos post my code, if needed.

You don't exactly need the whole Fast Fourier Transform (FFT) "mechanism" to achieve what you are trying to do. Also, your biggest problem in this setup is going to be synchronisation. Here is why:

You are looking at something like this:

Source --> Low Pass Filter (Extracting the bass) --> Rectification --> Integration --> Threshold --> Motor

You can implement this setup with discrete convolution, which means that you can have this pipeline operating on a per sample basis. A new sample comes in, it is processed and as a result a new processed sample pops out of that pipeline.

At its most basic and simplest form, a digital filter is a set of coefficients that are multiplied with the different time instances of the incoming signal in a rolling window fashion. Therefore, the "problem" now is how do you calculate those coefficients. There is a very large number of methods to achieve this but the simplest way to design a low pass filter is a sinc windowed filter. For more information, please see this link.

So, now you have your filter coefficients. Let's say you have designed a low pass filter with 50 coefficients and a cut off frequency of 100Hz. The filter is returning a signal that contains ONLY the first 100Hz of the Source.

How do you use this to drive the motor?

If you are not doing Pulse Width Modulation (PWM) on the motor, to create different intensities of vibration, then all that you are after is an "on-off" signal.

This "on-off" signal can be produced by a threhsold function. In other words, if the signal becomes "louder" than some pre-set value, then turn the motor on, otherwise, keep the motor turned off.

To judge if the low pass filtered signal is loud enough, you need to first rectify it (take its absolute value) and then integrate it (apply yet another low pass filter to smooth the data). This is equivalent to classical rectification and smoothing out. The low pass filter for smoothing here can be a sinc filter as above, with a very low cut off frequency (so that it becomes an integrator). Alternatively, you could use a moving average filter. For a given length, say $N$, the moving average filter's coefficients are $h_i = \frac{1}{N}, i \in [0,1, \ldots, N-1] $.

At the end of this process, you will have a "pulse" that goes high when the region below 100Hz is "loud" and low otherwise.

Finally, apply a threshold to this pulse (for example $x>T$ where $T$ is some value in the range of the signal and you have your square pulse that drives the motor.

You now have to put this "on-off" signal in a tiny little data packet and send it via serial (or other way) to the Arduino which will be doing the hardware bit.

In fact, if you were to use fixed point arithmetic, you can handle part of this process entirely in the Arduino. Unfortunately, they don't go "audio quality" high in sampling so that you could carry out everything in the Arduino.

Your biggest problem here is going to be synchronisation. Due to the way that digital filters work, there is a delay between the input and the output. The longer your (FIR) filters are, the better filtering they will be doing but also, the longer the delay they will be introducing to your signal. On top of this, you are going to have to add the delay of the connection with the arduino and the "packet encoding and decoding" that has to take place anyway.

In other words, your "kick" will be occuring slightly later than the bass and this might be annoying, depending on your specifications. You would have a similar issue even if you were implementing this with analog filters.

If you want the most absolutely low delay between input and output then you are going to have to move towards more complex to design filters.

Hope this helps.

EDIT:

To drive the motor with PWM and produce different levels of vibration for different levels of "bass" in the signal, you need to substitute the stages of "recitication --> integration" with a simple "quantisation" stage.

Quantisation is simply a linear transformation from the range of possible values of your "bass" signal (say 0-255) to the range of possible PWM levels (say 0-8).

The fact that the signal can have practically infinite dynamic range (AS COMPARED TO THE VIBRATION LEVELS) does not mean that the PWM stage can produce infinite vibration levels with the differences between them being discriminable.

Just as it happens with joystick controllers, the levels of vibration are usually very small, for example 4 or 8. There is no point in driving the motor in one instance with 2.812721 ms pulses and in the next with 2.812728ms pulses. The difference is incredibly small to be felt.

Quantisation can be as simple as round(bassLevel/maxBassLevel*vibrationLevels) to produce a number between 0 and vibrationLevels which you can then map to PWM periods through a lookup table. For more information on quantisation, please see this link

So, now your pipeline is:
Source --> Low Pass Filter (Extracting the bass) --> Quantisation --> Motor

And this will be relatively faster in responding because you got rid of one low pass filter.


How to Reduce Speech Feature Dimensions? I have been experimenting with speaker independent automatic speech recognition.  At this point, I understand that the signal is typically segmented into overlapping frames and MFCCs are extracted as features for each frame.

A lot of the texts I read talk about reducing the dimensionality of the features using (for example) Linear Discriminant Analysis or Primary Component Analysis.  While I think I understand the concepts involved in reducing dimensions in the abstract, I am not understanding how this feature reduction is supposed to be applied to the MFCC vectors.

I am interested primarily because several texts seem to suggest that reducing the dimensionality produces less variant features, which make recognizing patterns easier.  As a bonus, there is also less computational overhead.

Is the dimension reduction supposed to be applied on a per-frame basis?  For something like LDA, don't you need to know what class the data belongs to?

Thanks,
Carl

You have one feature vector per frame. The collection of feature vectors for all the frames in the training data form a matrix. Since this is training data, you know the true class for each frame. Use this matrix and the training data labels to build your LDA matrix. Alternatively, just use this matrix (without the labels) to find the PCA matrix.

Now you use this matrix to reduce the dimensionality of the feature vector for each frame in your test data when you're doing the speech recognition.

PCA doesn't produce less variant features. I think what you meant to say was that the features are less correlated with each other. This is called whitening and it can help some classifiers better learn from the training data.


How do KLT and SVD differ from each other? What are the differences between KLT and SVD? It looks more or less similar to me. Did I misunderstand the concepts? Explain me more about SVD and if possible provide me the link for SVD (Image Processing). Thank you.

From the direct implementations of the two, SVD and KLT are different since $S$ in SVD is a diagonal matrix, while $C$ in KLT is not. Besides, KLT is originally used to describe a stationary process. In its expansion, an infinite number of basis functions is required to
form a complete set. Yet for a random vector $X$, it can be rephreased in terms of a linear combination of orthonornal basis vectors, then covariance matrix $C$ is furtherly derived.

For a 2D image, we can view it as $M$ realizations of stochastic process $x$ where $M$ is the row of the image. And the assumption is that all rows have the same row covariance matrix $C$, then the KLT can be somehow converted to the form of SVD. In other words, when considering a multivariate stochastic process, the SVD and the KLT are computationally equivalent. That's why in most cases people view the both as the same. For mathematical details, refer to The KL transform and eigenimages.

In the application of image compression, the global KLT is used for the whole image, while SVD is better applied to each block. In Patrick et al's Hybrid KLT-SVD Image Compression, they combined the two method to implement the adaptive transform on the blocks to achieve a better coding efficiency.


Gaussian pyramid: why needs the image to be downsampled In the context of a gaussian pyramid, why is the image downsampled separately although the numbers of pixels are decreased through smoothing already?
In my understanding, the number of pixels will decrease due to the smoothing process already, if a gaussian filter is used within the image boundaries. An image with 8 x 5 pixels smoothed with a gaussian filter of 3 x 2 pixels will create a 6 x 4 pixel image if the filter can only be slided within the boundaries of the original image without crossing the edges.

So if the constraint of not crossing the edges of the image is applied, why would one use additional downsampling?

In the context of a gaussian pyramid, why is the image downsampled separately although the numbers of pixels are decreased through smoothing already?

The Gaussian filter adjusts the bandwidth of the content of the image. The downsampling adjusts the spatial resolution of the image.

The Gaussian filter is a low pass filter. After its application to an image, only the frequencies within its pass band will be present in the resulting image.

An image of dimensions $M \times N$ pixels can depict $\frac{M}{2} \times \frac{N}{2}$ lines per pixel of spatial frequencies.

If you apply a low pass filter at $\frac{M}{4}, \frac{N}{4}$, you have reduced the spatial frequency requirements for that signal by half. Where before you needed $\frac{M}{2} \times \frac{N}{2}$ lines per pixel, now you know that spatial frequencies in your image do not exceed $\frac{M}{4}, \frac{N}{4}$ lines per pixel. The rest of the bandwidth that an image of that resolution can support goes unused.

At that point, you can reduce the image resolution to accommodate the "new" limits of the bandwidth without any loss.

From another point of view, you use the low pass filter to ensure that the resulting image will not contain any spatial frequency that could be aliased if it was to be plainly subsampled (because ultimately, resizing an image means that it will have to be resampled and "fit" into a smaller space. This smaller space means less bandwidth). Once you ensure that no frequency will be aliased, then it is safe to resize.

  An image with 8 x 5 pixels smoothed with a gaussian filter of 3 x 2 pixels will create a 6 x 4 pixel image if the filter can only be slided within the boundaries of the original image without crossing the edges.

This is to do with choosing only the "valid" values in the output of the filter. Values for which all the weights of the two dimensional impulse response are applied to image pixels.

Hope this helps.


HMM-based speech recognition/synthesis: state representation of monophones and triphones In literature, a monophone is represented with 3 states.

When doing context dependent clustering, 3 states are used to represent a triphone.

Why is the triphone represented with 3 states and not 5? Doesn't modeling with 3 states reduce the information?

My current understanding is that each of the phones in a triphone corresponds to each of the states. For example, given a triphone 'k-a+b', then the first state corresponds to 'k', etc.

My current understanding is that each of the phones in a triphone corresponds to each of the states.

That's not correct.

Each phone, regardless of the context dependency, has three states: "start", "in-between", and "end". There are always transitions from the previous state to the next one and to itself. This defines the structure of HMM. Now context-dependency means that for each triphone we know its left and right context (context-independent phone). Let's take a look at what phone AX could look like in both CI and CD HMMs for the word "one [W AX N]".

When context-independent phones are used AX is begin split it into three states - _AX_B_, _AX_I, _AX_E - to construct the corresponding HMM:

In case of context-dependent phones AX actually doesn't exist any more. Instead we have many new phones that store information of phones to the left of it and to the right of it. Thus in the given word the neighbours of AX are W(left) and N(right), so HMM for phone W-AX-N is constructed:


Amplitude of the Signal in Frequency Domain different from Time Domain How does the Amplitude of the Signal is changed when taking the FFT of a Signal,

Have a look,

The amplitude had changed from 2 to 30,

and Here is my code for generating the above output,

f=1e3;   %Frequency of Wave
A=2;    %Amplitude
Fs = 1e6;  %Sampling Frequency

Ts = 1/Fs;   %Sampling Rate

t = 10/f;       %Time period of 1 Oscillation = 1/f

n = 0:Ts:t;  %Generating Samples

x=A*sin(2*pi*f*n);
subplot(2,1 ,1)
plot(x

% 100,000 = Fs
% 10,000 = Length of the Signal
% 100,000/10,000 = 100Hz <- First point in FFT = 100Hz
% 2nd Point = 200Hz
% 3rd Point = 300Hz
% 4th Point = 400Hz
% .
% .
% 10th Point = 1KHz <- Original Signal Frequency

subplot(2,1 ,2)
F=fft(x);
plot(real((F))),grid on
xlim ([0 20])

You're overlooking four things:

The $\frac{1}{FFT\_size}$ normalization coefficient. Some FFT implementations have or do not have this factor. Check the definition of FFT as performed by matlab on the Mathworks site!
Why are you looking at the real part only? The amplitude is conveyed by the modulus (magnitude) of the complex number. Here, the real part is 31 but the imaginary part is -10000.
A sine wave of amplitude 2 yields two complex exponentials of magnitude 1.
Keep in mind that you are not computing the Fourier transform of a sine wave (which goes from minus infinity to infinity), but of a sine waved multiplied by a rectangular window. As a result, you won't find the spectrum of a sine wave, but the spectrum of a sine wave convolved by the spectrum of your aperture (a sinc function for the rectangular window). The energy will bleed over adjacent frequency bins instead of being confined in bin 11. You will need to sum the energy over the adjacent bins to recover all the energy!

Getting back to your example:

abs(F(11)) is 10000. Divide by the sample length and you get a bit less than one (which is the expected value). Sum the small residues across the adjacent bins and you'll get back your expected energy.


Is there a FFT algorithm with the circular buffering? I have an audio signal and i want to process it sequentially. Is there an algorithm which allows to modify the spectrum of a signal when i remove a part of signal and add the other one?
For example, i obtain a set of 256 samples of input data. Then i do 256 FFT and i want to add the result spectrum to 16384 FFT which is the big buffer. So i have the  big FFT buffer with the good frequency resolution and i should sometimes modify to correspond the spectrum as if i do 16384 FFT.
The main objective is i need to avoid 16384 FFT because it's quite expensive. The way i want to optimize is (in time domain):

Shift the buffer to 256 samples;
Place the new samples to last 256 entries in buffer;
Make 16384 FFT.

So whenever the 256 sapmles come i need to recompute 16384 FFT whereas i've already have 16384 - 256 computed samples. Is there way perform such task (circular buffering FFT)?
My thoughs was (in frequency domain):

Multiply the 16384 FFT buffer by the complex sinusoid (that coresponds to circular shift the signal in time domain to 256 samples);
Calculate 256 FFT of new portion of signal;
Expand the spectrum to 16384 by inserting zeroes between bins (it's corresponds duplication of time domain signal);
Mix the signal with removing duplication of the expanded spectrum (by convolving with sinc function), and zero last 256 samples (by convolving with sinc function).

Is there way to avoid two comlementar convolving in the last step?

The listed algorithm seems to be form of sliding DFT, summing narrow rectangular slice windows.  In general, an FFT (O(NLogN) is faster than a sliding DFT of slice size b (O(b*N)), for b > LogN , due to a factorization that is not available to slices alone.  Thus, avoiding the 16384 FFT seems to be more expensive than just doing the FFT (if you actually want all 8192 correct bin results), since 256 > log(16384).


Defining an average in time for a density in time I'm reading the text book "Time-Frequency Analysis" by Leon Cohen and I've made my way through a decent portion. There is however a conceptual issue I keep coming back to. The book states:

  If we consider $\lvert s(t)\rvert^2$ as a density in time, the average
  time can be defined in the usual way any average is defined:

  $$  \langle t \rangle  = \int t \lvert s(t)\rvert^2dt$$

I don't understand conceptually why multiplying by $t$ in the integration will provide us with an average. I didn't know that was the typical way to define an average either despite the books insistence that this is common knowledge (which it very may well be I'm just not interpreting it correctly). I've always thought about an average as dividing by some base, say by the total period of time $T$ NOT multiplying by $t$ inside the integral. It states earlier that:

  $\lvert s(t)\rvert^2$ is the energy per unit time.

So I was thinking that perhaps mathematically multiplying by $t$ which has units of $s$ would get rid of the per unit time portion, but the units of $\lvert s(t)\rvert^2$ aren't actually $v^2/s$ as that would be the derivative, it's just $v^2$ so that doesn't make much sense (assuming the signal is an electrical signal and the unit time is seconds). Can someone point me in the right direction on how to think about this? I know it's a simple concept but I'm just having trouble making it over this hurdle. Thanks!

I think it is implied that $\int_{-\infty}^{\infty}|s(t)|^2 dt = 1$ (which is why the author refers to $|s(t)|^2$ as a density).

I agree with your sentiment "I've always thought about an average as dividing by some base", but in this case, the base would be 1, so no division is needed. If you were to plot $|s(t)|^2$ as a function of $t$, and think of this plot as a histogram, the the integral:

$$\int_{-\infty}^{\infty}t|s(t)|^2 dt $$
will tell you the average value, since each value of $t$ is weighted by the density function.


Plot real-time FFT using matplotlib I'm trying to plot fft in python. I use pyalsaaudio for capturing audio in PCM (S16_LE) format. I use the ion() and draw() functions in matplotlib to have the fft plotted in real time. This is the program I wrote :

import alsaaudio as alsa
import numpy as np
from matplotlib import pyplot as plot
from matplotlib import animation
import time

#Configuration

card = 'default'
audio = alsa.PCM(alsa.PCM_CAPTURE,alsa.PCM_NONBLOCK, card)

def configure():
    plot.ion()
    audio.setchannels(1)
    audio.setrate(44100)
    audio.setformat(alsa.PCM_FORMAT_S16_LE)
    audio.setperiodsize(1000)

def run():
    loops = 100000
    plot.show()
    while loops > 0 :
            loops-=1
            length,data = audio.read()

            if length:
                #converting into a float array
                npdata = np.fromstring(data,dtype='<i2')
                #Obtaining FFT
                freq_list = np.fft.fft(npdata)
                plot.plot(freq_list)
                plot.draw()
                plot.clf()

    plot.close()

configure()
run()

1) The entire plot is jumpy. The y axis limits range from -5 t 5 to several tens of thousands. Since the range of y axis keeps changing, so does the location of the plot. How can I fix this?

2) The x axis range is the argument to audio.setperiodsize(). In this case, it is 1000. From the pyalsaaudio documentation, "When the hardware processes data this is done in chunks of frames. The time interval between each processing (A/D or D/A conversion) is known as the period". Hence period is something to do with the time domain. It should have no relevance in the plot since fft is in the frequency domain. So why does the x axis of my plot end at the period size?

3) I believe that the fft function returns a set of complex numbers. What do they mean? The real part signifies the frequency and the imaginary part the amplitude? When I call matplotlib.plot() on the set of fft values, what is being plotted actually?

4) My objective is to determine the dominant frequency of a signal. To do this, I need to know the frequency and the amplitude of the component with that particular frequency. Can I get this from the fft plot?

Thanks for looking into this.

Since the range of y axis keeps changing, so does the location of the
  plot. How can I fix this?

Use ylim to set the range of values on the Y axis.

  So why does the x axis of my plot end at the period size?

Your input vector is of size 1000 ; so the FFT result is a sequence of 1000 complex numbers. Since you don't specify values / units for the x axis of your plot, you just get an arbitrary index.

To get a frequency scale, label your x axis as arange(0.0, length - 1) / length * 44100

  What do they mean? The real part signifies the frequency and the imaginary part the amplitude?

You should really take some time to understand the Fourier transform. The use of freq_list as a variable name makes me think that you are mistaken as to how the FFT works.

The complex number in freq_list[i] represents the complex amplitude at frequency i / length * 44100. The magnitude of a complex amplitude represents actual amplitude, the angle of a complex amplitude represents phase.

You're probably only interested in the (squared) magnitude, and you'll get much nicer plots if you plot it on a logarithmic scale.

  My objective is to determine the dominant frequency of a signal. To do this, I need to know the frequency and the amplitude of the component with that particular frequency. Can I get this from the fft plot?

Yes, by looking at the index of the FFT bin with the largest magnitude, but it's a rather terrible method for determining the fundamental frequency of audio signals - it has poor resolution and will provide incorrect results in a variety of common audio signals (for example it is common for speech signals to have lower energy at their fundamental frequency than at one harmonic).


What is the Fastest Way Ever to Do 2D Correlation? Here is what i gathered so far:

Any sliding window classification, image filtering or similar can be fastly done by a FFT (flip the signal and do convolution).
Also, if the template/filter kernel is separable it is possible to do fast correlation by simply separating into multiple kernels and applying then sequentialy.
I think it is possible to do fast correlation by using integral images, but I am not sure how.

The question is: is FFT the fastest way? If now, which is? Consider both GPU and CPU implementations

If it is not: is there algorithm capable of separating a non-separable kernel (in aproximated way, logically)?

Context:

My filter is an Inner Product Detector, which is found by

$$
h = R^{-1}\mu_1
$$
where $ R^{-1} $ is the autocorrelation matrix and $\mu_1$ is the mean of a given class.

The matrix $ R $ is hermitian. And can be calculated by

$$
\sum x^Tx
$$

where $x$ are d-dimensional vectors and $[.]^T$ denotes the transpose.

Also: kernels are squares of size 21,23,25,27 or 35
Images are 128x128

In your case, since you have multiple images while you have a given set of kernels the DFT based Correlation would be the best fit.

Pay attention that the DFT Based Convolution / Correlation Is Equivalent to Convolution / Correlation with Periodic / Circular Boundary Conditions.
It means that if you need different boundary conditions (Like padding with Nearest Neighbor) then it means you need to pad the image large enough so the cyclic effect won't hurt you.

In that case you might be better doing, given enough memory, what's being don in DNN.
Apply the Correlation / Convolution with Matrix Multiplication.

Matrix Multiplication is so optimized in our days that even if its complexity is higher, in real world it is faster in many cases.
So you build the Convolution Matrix (You can even chain them) and the images matrix and multiply. I can show a MATLAB Example.

By the way, even if you do it with DFT in the GPU, I suggest you do it with a DNN library (Something which utilize CUDNN, ROCm - Chainer, PyTorch, TensorFlow, MATLAB Neural Network Toolbox, Knet, etc...) since this case of multiple images (Batch) with given kernel is the case optimized in those libraries with automatic mechanism to chose between DFT based or Matrix based implementation.

Good reference - Algorithms for Efficient Computation of Convolution by Karas Pavel and Svoboda David.


Why does the BER rate change randomly versus the SNR? I have altered a MATLAB code that simulates the BPSK modulation using a raised cosine filter at the transmitter and  a matched filter at the receiver to study how the signal changes in each step. But, the final plot of the BER versus SNR is not decreasing as the SNR increases. I have no clue how to fix this. I would appreciate any help. The image and code is below

BER vs SNR

clear
N  = 10^6; % number of bits or symbols
T  = 1; % symbol duration of 1us
os = 5; % oversampling factor
fs = 5/T; % sampling frequency in MHz
rolloff = 0.05;

Eb_N0_dB = [0:10]; % multiple Eb/N0 values

for ii = 1:length(Eb_N0_dB)
    filter = rcosine(1/T,os,'sqrt',rolloff);
   % Transmitter
   ip = rand(1,N)>0.5; % generating 0,1 with equal probability
   s = 2*ip-1; % BPSK modulation 0 -> -1; 1 -> 1

   % up sampling the signal for transmission
   sU = [s;zeros(os-1,length(s))];
   sU = sU(:).';
   sFilt = 1/sqrt(os)*conv(sU,filter);
   sFilt = sFilt(1:N*os);

   % Noise addition
   y= awgn(sFilt,ii,0);

   % mathched filter
   yFilt = conv(y,fliplr(filter)); % convolution
   ySamp = yFilt(os:os:N*os);  % sampling at time T

   % receiver - hard decision decoding
   ipHat = real(ySamp)>0;

   % counting the errors
   nErr(ii) = size(find([ip- ipHat]),2);

end

simBer = nErr/N; % simulated ber
semilogy(Eb_N0_dB,simBer,'mx-');
grid on
xlabel('Eb/No, dB');
ylabel('Bit Error Rate');
title('Bit error probability curve for BPSK modulation');

There are a few problems with your code:

You are using the index ii not the actual SNR in DB in this line: y= awgn(sFilt,ii,0);
You are not taking into account the filter delay in ySamp = yFilt(os:os:N*os);.  Your filter is filter = rcosine(1/T,os,'sqrt',rolloff); which will design a root-raised cosine FIR filter.  Your filter delay will be (length(filter)-1)/2.