Why the RMS of a PSD curve is the root of the area below I will try to explain what is my level of understanding of this problem, please correct me if I'm wrong:

RMS is the Root Mean Square, it represent the mean value of the input signal.
PSD is the measurement of the responses that shows me at which frequencies most of the energy is concentrated.
The area below a curve is the integration of that function.

My situation is that several random vibration tests are performed. These tests are called random tests because of the input signal. In contrast to a sine test where the structure is excited with a sinusoidal input, only one frequency is excited at a time, here 'all' frequencies are excited at the same time.

In this case PSD is measured in ${{g^2}/{Hz}}$ and RMS in ${g_{RMS}}$. Armed with that it easy to see that if you multiply PSD per the frequency range and you take the root of the result you will get something in ${g}$'s, but I don't know how exactly derive the famous relation:

\begin{equation}
{g_{RMS}=\sqrt{\int_{f_1}^{f_2}PSD(f)df}}
\end{equation}

The understanding that I have is very basic and it would be great if someone give me a clear idea of the relations among the RMS, PSD and the real signal. Thank you very much.

In the figure I have plotted a standard random signal.

I think this is simply an aspect of Parseval's Theorem (e.g. click me)

It simply says: sum of squares in the time time domain equals sum of squares in the frequency domain. Substitute "sum" for "integral" if using the analog domain. In other words: total energy in the time domain equals total energy in the frequency domain.

This can easily reproduce your formula, $g_{RMS}$ represents the  time domain energy and the integral on the right represents the energy in the frequency domain. The exact scaling depends on the details like length of the signal, periodicity, sample rate, etc.

I think you mean why is the rms of the voltage signal in the time domain the square root of the area below the PSD. This is because the area below the PSD is the average power and average power is the square of V RMS. The RMS of the PSD curve itself is in W/Hz and is not the square root of the area under it which would be in $\sqrt{W}$ i.e. $V$
Firstly, $X(2\pi)$ is the complex coefficient at the frequency of $2\pi$ i.e. the sum of [the area under curve of [all real and imaginary (co)sinusoids at the frequency $2\pi$ with their individual phases and amplitudes times a complex sinusoid at frequency $2\pi$]] in $Vs$
$$\frac{1}{2\pi}\int _{-\infty}^{\infty}|X(\omega)|\,d\omega$$
Is the (continuous) sum of the magnitudes of the complex coefficients of all frequency components of the signal, which results in the average complex magnitude divided by the time interval $[-\infty, \infty]$ in $V$. This is the average voltage of [all real and imaginary (co)sinusoids in the signal with their individual phases and amplitudes and frequencies times a complex sinusoid at those frequencies]
Squaring this you get a power value but it's a meaningless one.
The sum of the average voltage of the product of each real and imaginary (co)sinusoid in the signal with their individual phases and amplitudes at a specific frequency times a complex sinusoid at that frequency is:
$$\lim_{T\to\infty}\frac{1}{2T}|X_T(\omega)|$$
The ASD is this sum of the average voltage at a specific frequency divided by the square root of the bandwidth:
$$A_{xx}(\omega) = \lim_{T\to\infty}\frac{1}{\sqrt{2T}}|X_T(\omega)|$$
Anyway,
$$E = \frac{1}{2\pi}\int _{-\infty}^{\infty}|X(\omega)|^2\,d\omega \equiv \int_{-\infty}^{\infty} |x(t)|^2~dt$$
Is the (continuous sum) of the absolute square of the complex coefficients of the components of the signal. I.e. you get the average of the absolute square of the complex coefficients divided by the time interval $[-\infty, \infty]$ in $Ws$ or $J$. It is the energy of the signal.
$$P_{xx}(\omega) = \lim_{T\to\infty} \frac{1}{2T}|X_T(\omega)|^2$$
Is the absolute square of the magnitude of the total complex coefficient at that particular frequency divided by the time interval $[-\infty, \infty]$. It is the PSD, like $E$ it is also in $Ws$ or more appropriately $W/Hz$. The difference with $E$ is that $E$ is the average $Ws^2$ (ESD) divided by $T$($s$) whereas PSD is the $Ws^2$ at a specific frequency divided by $T$($s$). The area under the PSD is the power within a particular bandwidth.
$$\lim_{T\to\infty} \frac{1}{4T^2}|X_T(\omega)|^2$$
Is the sum of the square of the average voltage of the product of each real and imaginary (co)sinusoid in the signal with their individual phases and amplitudes at a specific frequency times a complex sinusoid at that frequency.
$$P_{\text{Avg}} = \frac{1}{2\pi}\int _{-\infty}^{\infty} \lim_{T\to\infty} \frac{1}{2T}|X_T(\omega)|^2 \,d\omega \equiv \lim _{T\to \infty }{\frac {1}{2T}}\int _{-T}^{T}|x(t)|^{2}\,dt$$
Is the continuous sum of [the absolute square of the magnitude of the complex coefficient at that particular frequency divided by the time interval $[-\infty, \infty]$] i.e. it is the average energy of the complex coefficients of the signal divided by the time interval $[-\infty, \infty]$, which is the same as the average power of the signal.
Take $5+\cos(\omega t)$. The average power of the signal in the time domain is $25.5W$ (the average of $(5+\cos(\omega t))^2$). The frequency domain transform is $π (δ(ω - 1) + 10 δ(ω) + δ(ω + 1))$. $|X(2\pi)|^2$ has the units $V^2s^2$. One of the frequency domain impulses has a magnitude of $V\times 2T$. When this is squared, it can be represented as an impulse with a magnitude of $V \times 2T \times V \times 2T$. An impulse of this magnitude would have an area of ${V}^2\times 2T$. When this area, which expands with $T$ in either direction is divided by $2T$, and the limit to infinity of this relationship is taken, you get the $V^2$, which is $\frac{2\pi(0.5^2 + 5^2 + 0.5^2)}{2\pi} = 25.5W$.
The complex coefficient itself is the average real and the average imaginary voltage of [all imaginary and real (co)sinusoids in the signal at that frequency (which have their own individual phases and amplitudes) times the complex sinusoid at that frequency] times the integration interval.
The continuous-time periodogram is (where $2T$ is the length of the symbol):
$$S_{xx}(\omega) = \frac{1}{2T}|X_T(\omega)|^2$$
Which is used for continuous aperiodic signals that have a certain pulse width where the rest of the signal to infinity is 0, like an incoming OFDM symbol while it is in the analogue domain. The interval of integration in this case is finite and the frequency domain presents sincs at those frequencies instead: $A\tau\frac{\sin\left(f\pi \tau\right)}{f\pi\tau}$, where $A\tau$ on the real sinc is the sum of the real parts of [the areas in $Vs$ under all the real and imaginary cosinusoids and sinusoids in the pulse window at that frequency times a complex sinusoid of that frequency] and the $A\tau$ on the imaginary sinc is the
sum of the imaginary parts.
The discrete-time finite sample periodogram is typically:
$$S_{xx}(\omega) = \frac{1}{N}\left|\sum_{n=0}^{N-1} x_n e^{-i\omega n \Delta t}\right|^2$$
This can be reformulated as:
$$S_{xx}(\omega) = \frac{1}{T}\left|\sum_{n=0}^{N-1} x_n e^{-i\omega n \Delta t} \Delta t\right|^2$$
Which is a sum of the real parts (and a separate sum of the imaginary parts) of [a riemann sum approximation of the area under the envelope created by [the samples of each real and imaginary (co)sinusoid in the signal at the frequency inserted into the function (with their own amplitudes and phases) times the complex sinusoid at that frequency]], and the output is a complex number with a real part of the result of the real sum and an imaginary part of the result of the imaginary sum, whose complex modulus is then squared and divided by the number of samples.
The power of a finite sample discrete-time signal is
$$P_{\text{Avg}}= \frac{1}{N}\sum_{m=0}^{N-1} \hat{S}_{xx}(m\Delta \omega) =\frac{1}{N}\sum_{m=0}^{N-1} \frac{1}{N}\left|\sum_{n=0}^{N-1} x_n e^{-im\Delta\omega n \Delta t}\right|^2$$
$$\equiv \frac{1}{N^2}\sum_{m=0}^{N-1} \left|X[m]\right|^2 \equiv \frac{1}{N}\sum_{n=0}^{N-1} \left|x[n]\right|^2 $$

For a wide-sense stationary complex random process:
$$P_{xx}(\omega) = \lim_{T\to\infty} \frac{1}{2T}\mathbf{E}\left[|X_T(\omega)|^2\right] = \lim_{T\to\infty} \frac{1}{2T} \mathbf{E} \left[\int_{-T}^{T} x^*(t) e^{i\omega t}\, dt \int_{-T}^{T} x(t') e^{-i\omega t'}\, dt' \right] = \lim_{T\to\infty} \frac{1}{2T} \int_{-T}^{T} \int_{-T}^{T} \mathbf{E}\left[x^*(t) x(t')\right] e^{i\omega (t-t')}\, dt\, dt'.$$
So we can transform this into:
$$\lim_{T\to\infty} \frac{1}{2T} \int_{-T}^{T} \int_{-T}^{T} \mathbf{E}\left[x^*(p) x(p+\tau)\right] e^{-i\omega \tau}\, dp\, d\tau.$$
i.e.
$$\lim_{T\to\infty} \frac{1}{2T} \int_{-T}^{T} \int_{-T}^{T} R_{xx}(\tau) e^{-i\omega \tau}\, dp\, d\tau.$$
By substituting $\tau = t'-t$
and $p = t$ and the determinant of the jacobian matrix that transforms $dt~dt'$ to $dp~d\tau$ is $1$
The expected value removes the dependency on $p$ because it's a constant:
$$\lim_{T\to\infty} \frac{1}{2T} \int_{-T}^{T} \int_{-T}^{T} R_{xx}(\tau) e^{-i\omega \tau}\, dp\, d\tau.$$
$$=\lim_{T\to\infty} \frac{1}{2T} \int_{-T}^{T} 2TR_{xx}(\tau) e^{-i\omega \tau}\, d\tau.$$
$$=\lim_{T\to\infty} \int_{-T}^{T} R_{xx}(\tau) e^{-i\omega \tau}\, d\tau.$$
$$=\int_{-\infty}^\infty R_{xx}(\tau) e^{-i\omega\tau}\,d \tau=R_{xx}(\omega)$$
The continuous time periodogram is (where $2T$ is the length of the symbol):
$$S_{xx}(\omega)=\mathbf{E}\left[\frac{|X_T(\omega)|^2}{2T}\right] = \int_{-T}^{T} R_{xx}(\tau) e^{-i\omega \tau}\, d\tau.$$
Where we need the PDF to compute the ensemble average $\mathbf{E}$
The discrete time periodogram is:
$$S_{xx}(\omega) = \mathbf{E}\left[\frac{1}{N}\left|\sum_{n=0}^{N-1} x_n e^{-i\omega n \Delta t}\right|^2\right] = \sum_{\tau=0}^{N-1} R_{xx}(\tau) e^{-i\omega \tau}$$
Where we need the PDF to compute the ensemble average $\mathbf{E}$ but if $N$ is ergodic then a time average can be used. You would get:
$$\lim_{N\to\infty} \frac{1}{N-\tau}\sum_{\tau=0}^{N-1}\sum_{n=0}^{N-1-\tau} x[n]x^*[n+\tau] e^{-i\omega \tau}$$
The estimator of the periodogram of a WSS uses one record:
$$\hat{S}_{xx}(\omega) = \frac{1}{N}\left|\sum_{n=0}^{N-1} x_n e^{-i\omega n \Delta t}\right|^2$$


Why is this MATLAB butterworth filter giving me wrong results Given a 128 Hz signal, in MATLAB I entered the following command

[B,A] = butter(2,[1/64,63/64],'bandpass'); %construct second order all pass filter
Filtered_Data = filter(B,A, Data) %filter data using butterworth filter

And got

This is far from what I expected since the filter should have passed all frequencies. an someone help point out what I'm doing incorrectly?

Your filter is passing only between 1 and 63 Hz, assuming that you meant to say that the sampling rate is 128Hz.

Since you are removing everything below 1Hz, the plot you have looks reasonable.

The range of values in your original signal is relatively small (4395 to 4420,) and the range of values in the output plot is also relatively small (cannot tell exactly due to scaling, but it looks small.)

The high values (+4000, -2000) in your output plot are the step response of your filter.  Your data jumps from zero (before the first data point comes in) to around 4400 at the beginning.  This step causes the filter to overreact - the short peaks up to 4000 and down to -2000 are the result.  This damps out right away, and the rest looks pretty much like what would be expected.

The offset (about 4400) looks like a DC component to the filter.  Since DC is lower than 1Hz, it gets removed as well.

The problem might be in the way you are defining the cut-off frequencies.
Matlab uses normalized frequencies (look at the examples here http://it.mathworks.com/help/signal/ref/butter.html), therefore you have to define them in this way:

assuming you want to cut your signal at $5 Hz$ and $100 Hz$ and your sampling frequency is $Fs$ (which you haven't specify)

cut_f1 = 5;
cut_f2 = 100;
Wn = [cut_f1/(Fs/2) cut_f2/(Fs/2)];

This is the $Wn$ you have to pass to the butter() function.

Moreover, I suggest you to use the command

filtfilt()

to filter your data in such a way to null the phase shift introduced by the filter (http://it.mathworks.com/help/signal/ref/filtfilt.html).


Transmission Line Information Loss Suppose we have a communication channel that has a bit-error rate of 50% and the bits are flipped at random such that the receiver has no efficient way to tell which bits have been flipped and there are no error correction bits transmitted with the signal. How do you calculate the information loss (in bits) imposed by the channel? I know there a Shannon formula for doing this but I don't know what it is! Could someone solve this for me with a specific example? Thank you.

A binary input - binary output channel which flips bits with equal probabilities is known as a Binary Symmetric Channel.

(Image Credit: Wikipedia - Binary Symmetric Channel)

The channel capacity $C$ is the maximum (actually, a supremum) mutual information between the channel inputs and the channel outputs over input distributions. You can't get information through a channel more than its capacity. Intuitively, this says that by designing the way you put inputs into the channel well, you can maximize the amount of stuff coming out of the channel that you can recover (and this limit is the capacity).

Mathematically, this is
$$C = \sup_{p(x)} I(X;Y)$$
where $p(x)$ is probability distributions on the input $X$ and $Y$ is the output of the channel.

In this case, you can write the input of the channel as $P(X=0) = q$ and $P(X=1) = 1-q$. With this distribution, you can see that $P(Y=0) = P(X=0 \text{ and no crossover}) + P(X = 1 \text { and crossover}) = q(1-p) + (1-q)p$ and $P(Y=1) = 1-P(Y=0)$. The mutual information between $X$ and $Y$ is given by $I(X;Y)=\sum_{x,y} p(x,y) \log \frac{p(x,y)}{p(x)p(y)}$ where $p(x,y)$ is the joint distribution of $x$ and $y$, and $p(x)$, $p(y)$ are the distributions of $x$ and $y$ respectively. By writing down $I(X;Y)$ and then maximizing over $q$, you will see $C = 1 - H(p)$ where $H(p)$ is the binary entropy function $-p \log_2 p - (1-p) \log_2 (1-p)$. Then, in your case, $p = \frac{1}{2}$ where $H(p) = 1$, so the capacity is $0$.

If you're interested in learning more about this, a few good books to look (at very different levels of difficulty):

Elements of Information Theory (Second Edition), by Tom Cover and Joy Thomas (Easy, Probably Best Starting Point)
Information Theory and Reliable Communication, by Robert Gallegar (A bit harder)
Information Theory: Coding Theorems for Discrete Memoryless Channels by Czisar and Korner (Pretty hard)
Information Theory by Robert B. Ash (Dover book, Easy)

And an unrelated book which I love, but couldn't keep myself from mentioning, is David Mackay's Information Theory, Inference and Learning Algorithms (free on his website), even though it doesn't go over channel capacity.

the information capacity through a channel is, from Shannon,

$$ I = \int_0^B \log_2 \left( 1 + \frac{S(f)}{N(f)} \right) df $$

where $I$ is the rate of information passed through the channel in  bits per unit time, $B$ is the one-sided bandwidth, $S(f)$ is the power spectrum of the signal and $N(f)$ is the power spectrum of the noise.


Finding audio signal time delay: cross-correlation vs. impulse response? I need to find a time delay between the two audio signals, one comming from one set of speakers, and the other being the source signal that is initialy sent to speakers. I send some noise through those speakers and they come to microphone with a time delay compared to a source signal.

Cross-correlation method
As I understand, there is a cross-correlation method where I send both recorded signals to frequency domain via Discrete Fourier Transform (DFT), then multiply the frequencies, and finaly go back to time domain via Inverse DFT. There should be a peak in time domain telling me what the delay is (let's disregard the accuracy for now.)

IR method
There is also an Impulse Response method, where I also take two signals to freq domain via DFT, then divide the two to get a system transfer function, and finaly do an Inverse DFT of the Transfer function; this will give me the Impulse Response, which should also have a strong peak at the delay point.

My question is:
Am I right in the short descriptions above? The main difference between the two methods seems to be that frequencies are multiplied in the cross-correletion method, and dividied in the IR method? Have I missed something important?

Your description of the cross-correlation method is not entirely correct. What you describe is actually a convolution. For a cross-correlation you would perform an FFT on both signals and then multiply one with the complex conjugate of the other, and finally perform an IFFT.


Is there any Spread Spectrum system that uses Rake with cyclic Prefix? My question is that could the use of cyclic prefix can enhance the performance of Rake receiver ?

I find this proposal here CP for DSSS Rake Receiver that seems correct.

Without CP, ISI emerges and stays as interference in Rake finger.

With CP, no ISI.

The cost could be an overhead of transmission. But if I take a typical value of delay spread 5µs (ETU channel model) with UMTS SF = 256, the symbol duration is roughly SF/3.84Mcps = 256 / 3.84e6 = 66,66 µs. The cost 5/66.66 = 7.5% is acceptable ?

If the reason is not the overhead cost, is there any reason ? Or Rake receiver performance is not impacted by adding CP thanks to its pilots ?

you can read the following paper, they are using almost the same idea, the title of the paper: Symbol Cyclic Shift Equalization PAM-OFDM - A Low complexity CP-free OFDM Scheme ..  Good luck


Conditions for precoding matrix to preserve complex conjugate symmetry on DFT vector Suppose there is a DFT vector $\mathbf{X}$ with length N, which presents complex conjugate symmetry around its middle point, i.e., $X(1) = X(N-1)^*$, $X(2) = X(N - 2)^*$ and so forth.  $X(0)$ and $X(N/2)$ are the DC and Nyquist frequency respectively, therefore are real numbers. The remaining elements are complex.

Now, suppose there is a matrix $\mathbf{T}$, with size $N \times N$, which multiplies vector X.

\begin{align}
\mathbf{Y} = \mathbf{T}\mathbf{X}
\end{align}

The question is:

In what conditions, for matrix $\mathbf{T}$, the complex conjugte symmetry around the middle point of the resulting vector $\mathbf{Y}$ is preserved?

The motivation for this question is trying to come up with a precoder matrix $\mathbf{T}$ that results in a precoded (pre-equalized) symbol $\mathbf{Y}$ whose IFFT is real.

EDIT:

Thanks @MattL. and @niaren.
The difficulty about this question is to find necessary conditions. Matt's answer is indeed sufficient. It is also sufficient to make the following modifications:

The first row and the first column do not need to be zero. Instead, they could be non-zero, as long as its values present a complex conjugate symmetry around the middle point, its first value is real and its $(N/2+1)$-th value is real, just like the symbol. The same could be stated for the $(N/2+1)$-th column, the $(N/2+1)$-th row, and the main diagonal.

Secondly, the same correspondence between the matrix in the upper left corner and lower right corner could be made between the upper right corner and the lower left corner, that is, chose an $(N/2 -1)\times(N/2-1)$ matrix starting from $t_{2,N/2 + 2}$ to $t_{N/2,N}$, flip from left to right, flip upside down and take the conjugate, then put at the lower left corner. On MATLAB, that would be:

T(N/2+2:N,2:N/2) = conj(fliplr(flipud(Tisi(2:(N/2),N/2+2:N))))

This structure is similar to the structure of the DFT matrix. Would that be a necessary condition?

EDIT(2):

The following code implements such a valid operator for any real-valued $N \times N$ matrix $\mathbf{A}$:

N = 8;
A = rand(N,N); %must be real-valued
w = exp(-1j*2*pi/N); % twiddle factor
W = w.^(repmat(0:N-1,N,1).*repmat(0:N-1,N,1).'); % DFT matrix
T = W*A*W'

EDIT(3):

It is also interesting to note the $\mathbf{T}^{-1}$ presents the sufficient condition too. This comes from the fact that:

\begin{align}
\mathbf{T}^{-1} &= \mathbf{\left(W A W^{H}\right)}^{-1}\\
 &= \mathbf{\left(W^{H}\right)}^{-1} \mathbf{A}^{-1} \mathbf{W}^{-1}
\end{align}
where $\mathbf{W}$ is the DFT matrix.

Since $\mathbf{W^{H}} = N\mathbf{W}^{-1}$. This equation becomes:

\begin{align}
\mathbf{T}^{-1} &= \left(N\mathbf{W}^{-1}\right)^{-1} \mathbf{A}^{-1} \frac{1}{N}\mathbf{W^{H}}\\
&= \mathbf{W} \mathbf{A}^{-1} \mathbf{W^{H}}
\end{align}

Finally, since $\mathbf{A}^{-1}$ is real-valued, provided that $\mathbf{A}$ is full rank,$\mathbf{T}^{-1}$ is sufficient.

I think that the entries in your matrix $\bf{T}$ must obey $a_{N-n+1,N-m+1} = a^*_{n,m}$.
This is saying that the entries in row $N-n+1$ are the same as the coefficients in row n but where the coefficients are conjugated and reversed. The pattern in $\bf{T}$ for $N=4$ is

$T_4 = \left[ \begin{smallmatrix} a_{11}&a_{12}&a_{13}&a_{14}\\ a_{21}&a_{22}&a_{23}&a_{24} \\ a^*_{24}&a^*_{23}&a^*_{22}&a^*_{21} \\ a^*_{14}&a^*_{13}&a^*_{12}&a^*_{11}\end{smallmatrix} \right]$

I'm sure someone will come up with a better and more precise answer.

If I'm not mistaken the only solution for $\mathbf{T}$ which is independent of the vector $\mathbf{X}$ is a diagonal (complex) matrix, where the diagonal satisfies the complex conjugate symmetry.

EDIT: OK, I was mistaken. Diagonal is fine, but it's not necessary. The matrix $\mathbf{T}$ must have the following general structure: elements $t_{11}$ and $t_{N/2+1,N/2+1}$ must be real-valued (they correspond to DC and Nyquist). Apart from $t_{11}$ the first row and column contains only zeros. For elements $t_{22}$ to $t_{N/2,N/2}$ chose an arbitray $(N/2-1)\times (N/2-1)$ matrix. Then use this arbitray matrix to form a new matrix by swapping all rows (the first row becomes the last one, the second row becomes the second last one, etc.), by flipping the rows from left to right, and by conjugating. Then put this submatrix in the lower right corner of the total matrix $\mathbf{T}$. All other elements of $\mathbf{T}$ must be zero. I'm aware that this is a bit tough to understand without a visualization, so I'll add one later on when I have more time.


Choosing the best measure for machine learning on a Spectrogram I am working on an application to classify Talas(a piece of music with onsets). For that, I have the spectrogram of a Tala, which I want to compare with the test Tala, or compare 2 Talas.

For this I intend to find the onsets in the spectrogram and then form a classifier based on the frequencies that I get at the onset points in the Tala.

Now as we already know that a spectrogram contains multiple frequencies for a time t , I wish to quantify this as a single value to be used in machine learning, so that I have n values for the n onsets(say), and not multiple frequency values for the same time t.

So here is my question - Given all the frequencies at time t, how can I get a value that represents all these frequencies(like mean, median, etc). What should be the measure I should go ahead with?

First thing: a spectrogram does not give the "frequencies" at time $t$. It gives amplitudes, measured for a set of linearly spaced frequency bands.

Anyway, there is no "right" answer to your question, because it depends on which musicological aspect you want to focus on, and you need to enlighten us on that. If you want to focus on melody, extracting the fundamental frequency is the most relevant feature, but I suspect this is not the case here. If you want to discriminate different timbres, spectral centroid might be a better way of "summarizing" a slice of spectrogram into a single value - better, but still quite bad.

Indeed I can see several flaws in your methodology:

Why extract information at time t only (assumedly at a small window centered at t) rather than on the entire segment of sound between onsets?
Why do you have this constraint that for any onset you have to extract a single real value? Some properties of sound, like timbre, are multidimensional, so you won't be able to entirely capture them in one single value. If we take the example of discriminating tabla sounds (which might or might not relevant to your problem, since you mentioned tala...), we need several features to discriminate them (bayan/dayan discrimination relies on pitch information, while resonant vs non-resonant stroke discrimination relies on timbral features).
You are completely discarding the rhythmic (time) information. Is that really wanted for your problem?

Tell us more about your problem, there's a good chance it has already been solved, and that you are asking the wrong questions to solve it...

When working with audio signals there are two tricks, first remove the DC from the spectrum and band-limit it to a reasonable size. If you are looking at 1000 point intervals only use a 30 point FFT.

You can try to use ML to get the mean and variance and then use these vectors to evaluate test data. More advanced students might consider the MMI criteria.


displaying speech data in log scale using Matlab I have speech data that I can represent using matlab plot.
however, If i want to display it on log scale(like cooledit or Audition), I can take something like

data = audioread(pth);
for j=1:length(data)
    if data(j) >0
         data(j) = 20*log10(data(j));
    else
         data(j) = 20*log10(abs(data (j)));
    end
end

problem with this is that zero values map to infinity and negative values map to positive side . I need representation where Y axis would run in db like [0 -3 ..-21..-inf.. -21 -3 0](like in Cooledit) and the data is plotted as per this axis.
How do i do this?
sedy

It doesn't really make sense to display an instantaneous AC waveform on a log scale - as demonstrated by your results. I googled for cooledit screenshots and found this.

To accomplish this you would want to keep your data linearly scaled and draw it normally. Then just label the axis at the desired dB points.


Detect repetitive units within signals I have several signals that consist of repetitive units. In the figure you'll clearly see the variability of the signals, that increases top down. The first signal is super repetitive and units are indicated with green lines.
In the third, you'll see peaks and in the middle a little insertion that I know consists of rather diverging units, which still are units however (red lines).
The remaining three signals display the variation more.

Which signal processing / machine learning tool should I use in order to detect these units? With thresholding it works to find the significant peaks, but once signals get funky it's really difficult to accurately detect unit positions.

Edit:

I made some progress with spectral analysis and plotted a filtered signal over the initial signals. Indicated by the blue arrow is a region variable repeats. Here, the amplitude of the original and filtered signal do not match. Same for the red arrow. Similar effects for signal further down.

Your repetitive signal seem to be in a lower rate than the noise. I would run an FFT (why short) on all 6 sequences and filter out (zero) the frequencies that are clearly attributed to noise and then run an inverse FFT.
If you know the repetition rate you might select the one frequency that has the correct information and use the phase information to locate the desired pick. I will try to include exactly 20 cycles if possible (or 16 or 32).

ok, what I would do is, first take an autocorrelation of the give signal. you get the peaks(due to similar signal element overlap), find the time period by seeing the time scale difference of similar peaks.use the corresponding frequency to filter out the different signals.
does that help?


When negative frequencies cannot exist in practice, why do we calculate power in negative frequency? Suppose a signal with two-sided power spectral density $N/2$ is passed through an ideal LPF having cutoff frequncy  $B\,\mathrm{Hz}$. The output power is equal to $N/2 \times 2B = NB$ Watts. Why not $NB/2$?

Negative frequencies do exist. The vector base for Fourier transforms (and therefore power spectral density) is the set of complex exponentials $e^{j\omega t}$ where the frequency $\omega$ can be positive or negative.

A real sinusoid can be decomposed into complex exponentials:

$$\cos(\omega t) = \frac{e^{j\omega t} + e^{-j \omega t}}{2}$$

There you have your practical signal (the cosine wave) expressed in terms of the Fourier vector base (complex exponentials). Both sides of the spectrum must be considered in order to calculate power.

This is similar for any other real signal. In general, for real signals, the power spectrum will be symmetric and you must consider power in both sides of the spectrum.

The negative is comparatively concept .. this means when your reference signal has a main component @ fo there are also other components that leading and lagging to fo .. so the leading and lagging produce the (+) and (-) components on the right and left sides if fo respectively ..

I hope this may clear your confusion about negative frequencies..


Non Uniformity in MFCC Training Data I'm currently working on Person Identification, and using MFCC as the feature to my classifier. I have a database full of tagged(with labels) voices, and I use a CSV file to access the respective voice files. I generate the coefficients by dividing each voice file into equal frames of 1024 bytes & then, I find the coefficients of each frame. My problem is that all the training data isn't of the same size. Thus, the total number of coefficients vary from file to file. This dramatically affects the performance of my system. Is there any way around it? I've tried treating the whole file as a single frame, but the efficiency hasn't been very good. Any help would be appreciated.

Definitely don't use the whole file as a single frame.
As Alexander says, you want to use 25 millisecond frames with a 10 second overlap/hop (the vocal tract is stationary over this duration, so this is a common choice in speech recognition).
Each frame of audio makes a certain number of MFCCs (say 20). This number is constant for each frame of audio. But a file of different duration will, of course, have a different number of frames (and thus MFCCs). The size of your features (20) should always be uniform, but the number of observations (how many frames) you get from each file usually will not be, and that is okay.
Once you have extracted these acoustic-features, you need a machine learning algorithm that can "learn" what MFCC features are characteristic to a particular speaker. This is a very broad topic, but I can point to you this project: https://github.com/ppwwyyxx/speaker-recognition where you can see some of the machine learning algorithms used.

Good luck!

MFCC computation is based on spectral analysis, so you cannot feed an arbitrary signal to it in whole. The signal must be stationary. It's a common practice to assume that the speech signal is stationary within 25ms frame. Also it's a common practice to overlap frames by 10ms to detect transient events. Next, you cannot compare those MFCC directly. You rather should choose a pattern recognition method, such as DTW or HMM, where MFCC act as features.


Why images need to be padded before filtering in frequency domain In image processing books , we are told that Images need to be padded while doing filtering in the frequency domain. Why we need that zero padding ?

Multiplication in the frequency domain corresponds with circular convolution in the spatial domain.  This means that without padding the image properly, results from one side of the image will wrap around to the other side of the image.

You can think of 2D filtering as a sliding window that is centered over each pixel in the image and the center output pixel is a weighted sum of the pixels in the window.  With circular convolution, when the window hangs over the right edge of the image, it is actually wrapping back around to the left side of the image.  This means that output pixels on the right edge of the image will be affected by pixels on the left edge, which is almost never what is actually desired.

Zero-padding allows space for this wrap-around to occur without contaminating actual output pixels.

Fourier transform is periodic in nature.
Periodic function can cause interference between adjacent periods
and this will lead to wraparound error.

To overcome this we go for zero padding


Good book or reference to learn Kalman Filter I am totally new to the Kalman filter. I've had some basic courses on conditional probability and linear algebra. Can someone suggest a good book or any resource on the web which can help me can understand Kalman Filter operation?

Most websites start directly with the formula and what they mean, but I am more interested in its derivation, or if not detail derivation then at least the physical significance of each operation and parameter.

Many years ago I wrote this tutorial on the Kalman filter.  It derives the filter using both the conventional matrix approach as well as showing it's statistical assumptions as an 'optimal' least squares filter.

The best resource is Wikipedia page in my opinion. Here is a minimal and simple implementation of Kalman Filter with same notations given in Wikipedia page: https://github.com/zziz/kalman-filter


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

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

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

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

OK Lets go!

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

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

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

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

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

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

i might suggest looking up the good ol' Average Magnitude Difference Function (AMDF) or Average Squared Difference Function (ASDF) or Autocorrelation as means of determining the period which is the reciprocal of the fundamental frequency.


FIR filter design by the Fourier transform method I am having some problems understanding how the Fourier transform method is used to determine the FIR filter.

As far as i have understood, you start by using the ideal impulse reponse for the specific filter design you want to have.  Those are given in this table

As the the ideal impulse reponse contains infinite impulses, we truncate the infinite to finite a number by $m$, (which is confusing how this is done).

But my question is how do predetermine the order of the filter, without having to try randomly, using this method?

Some papers provide rule-of-thumb formula for the filter length (and coefficient quantization), for instance:

Accurate estimation of minimum filter length for optimum FIR digital filters, 2000, Ichige et al.
On computational complexity in digital filters, 1981, Bellanger
FIR Digital Filter Design Techniques Using Weighted Chebyshev Approximation, 1975, Rabiner et al.

The one I generally  use in lectures is borrowed from the second paper, for a low-pass filter with the following design:

Then, the estimated filter length $N$ (apparently $2M+1$ for you) is, with sampling frequency $f_s$:

$$ N = \frac{2 f_s}{3(f_2-f_1)}\log_{10} \left(\frac{1}{10\delta_1\delta_2}\right)$$

You design a filter with the desired frequency in mind that you whish to extract. So by designing such a filter you will elimnate undesired frequencies.


FFT of a AM modulated signal I am generating an AM modulated sinusoidal wave. The carrier frequency is set at $1000 \ \rm Hz$; the modulation frequency is set to $40 \ \rm Hz$ and $100 \%$ of the amplitude is modulated.
I generated 2 signals with 2 slightly different equations.

Signal 1 follows:
$$\sin(2\pi f_m t) \cdot \sin(2\pi f_c t)$$
while signal 2 follows:
$$\big(1-\sin(2\pi f_m t)\big) \cdot \sin(2\pi f_c t)$$

I don't understand why the first signal doesn't have a $1 \ \rm kHz$ component on the FFT and why both signals do not have the same period.
Both $\sin(x)$ and $1-\sin(x)$ have the same period, and I can see that signal 1 has the amplitude sinus and the carrier sinus in phase, but I can't get my head around those plots and interprete them correctly. Thanks for all the information and explanation you can provide :)

Code snippet to produce the 2 signals:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
import numpy as np
import os

#%% Signal 1
fs = 44100        # sampling rate, Hz, must be integer
duration = 1.0    # in seconds, may be float

# Define the time series
t = np.linspace(0, duration, int(duration*fs), endpoint=False) # time variable

# AM - Amplitude Modulation
fm = 40 # Modulation frequency
amplitude = np.sin(2*np.pi*fm*t)

# Output signal
fc = 1000 # Carrier frequency
signal1 = amplitude * np.sin(2*np.pi*fc*t).astype(np.float32)

# Apply FFT
fft_freq1 = np.fft.rfftfreq(signal1.shape[0], 1.0/44100)
fft1 = np.abs(np.fft.rfft(signal1))

#%% Signal 2
fs = 44100        # sampling rate, Hz, must be integer
duration = 1.0    # in seconds, may be float

# Define the time series
t = np.linspace(0, duration, int(duration*fs), endpoint=False) # time variable

# AM - Amplitude Modulation
fm = 40 # Modulation frequency
amplitude = np.sin(2*np.pi*fm*t)

# Output signal
fc = 1000 # Carrier frequency
signal2 = (1-amplitude) * np.sin(2*np.pi*fc*t).astype(np.float32)

# Apply FFT
fft_freq2 = np.fft.rfftfreq(signal2.shape[0], 1.0/44100)
fft2 = np.abs(np.fft.rfft(signal2))

#%% Plot
f, ax = plt.subplots(2, 3, sharex=False)
ax[0, 0].plot(t[:4411], signal1[:4411])
ax[0, 0].set_title('Signal 1')
ax[1, 0].plot(t[:4411], signal2[:4411])
ax[1, 0].set_title('Signal 2')

ax[0, 2].plot(fft_freq1[900:1101], fft1[900:1101])
ax[0, 2].set_title('Signal 1 FFT')
ax[1, 2].plot(fft_freq2[900:1101], fft2[900:1101])
ax[1, 2].set_title('Signal 2 FFT')

ax[0, 1].plot(t[:4411], amplitude[:4411])
ax[0, 1].set_title('Signal 1 AM')
ax[1, 1].plot(t[:4411], (1-amplitude)[:4411])
ax[1, 1].set_title('Signal 2 AM')

There two definitions for the AM modulated signals
The first one is called as the classical-AM (or conventional AM) and is given by
$$x_{AM}(t) = (A_c + m(t)) \cdot \cos(2\pi f_c t) \tag{1}$$
and the second one is called as the DSB-SC (double side-band suppressed carrier) and is given by:
$$x_{AM}(t) = m(t) \cdot \cos(2\pi f_c t) \tag{2}$$
In your code, you use Eq.2 (DSB-SC) for the first example, and Eq.1 (classical AM) for the second example. Their difference is that the first one includes the carrier spectrum (the middle spike at $f_c$) at the output in addition to the modulating signal spectra (side-band spikes). The latter DSB-SC, however, only includes the modulating signal side-band spectra
and not the carrier spectrum, hence the name suppressed carrier.
DSB-SC has the advantage that the tranmission energy is reduced as the permanently radiating carrier energy is saved yielding a more efficient transmission. Whereas the classical AM, which is less energy efficient to broadcast, has the advantage that its demodulation (by analog hardware) is extremely simple requiring what's known as the envelope detector (a diode, a capacitor and a resistance, in addition ot the antenna receiver and audio amplifier circuitry).

They do have a 1Khz component, but not via a Fourier lens.
To see what's happening, let's superimpose the AM signals with pure carrier tones:

Left plot reveals the key difference between the two AM schemes: polarity. Right is positive-only, left is zero-centered and flips signs. The effect of this sign reversal is, the "wiggling" of the carrier prematurely flips polarity. The red line shows when AM changes (+) to (-); without this sign reversal, we have a would-be positive peak at sample 32 (just as with blue pure tone) - with it, it's now negative.
This "fools" the Fourier decomposition into thinking that it's "more cycles", as that's what it would take to negative-peak sooner under a fixed-amplitude construction. The physical processes are still very much $f=2, 20$.
So the spectra we see are the result of having to represent the signal as a sum of fixed-amplitude frequencies. This Fourier construction is neither always the most meaningful, nor the most useful; I provide more intuition in this answer, and more on why FFT can mislead here and here; ignore the downvotes.

Caveat: if you count the orange peaks you won't find 40 in either case; each bears its own explanation given the AM. It may not be meaningful in the physical sense to decouple AM from carrier globally - e.g. consider what happens at right plot toward edges or right at center, where the AM is zero; an oscillation "skips", changing the "mean frequency" over the fundamental period. The left case shows more peaks, again per polarity flip. You can find more details here.
Further, you'll find slightly more cycles for orange even in the "well-behaved" region (away from AM zeros); it then may vary by purpose, but in an amplitude-decoupled framework this reflects as slightly higher carrier frequency, again more sensible than Fourier's. Lastly, arguably neither handle the case $f_m \rightarrow f_c$ too well in discerning what counts as "carrier" vs "modulator".

On sound: for predicting how we'd hear an AM waveform, a Fourier decomposition is indeed appropriate, as auditory processing is accurately modeled by fixed-amplitude sinusoidal decomposition (though not always), and the amplitude-decoupled decomposition won't apply.

Re: period there may be differences in fundamental periods $T$, but the two will still share a non-fundamental period. Put differently, the expected periodicity doesn't necessarily disappear, rather a smaller period is introduced due to AM to one of the signals and not the other. This is understood again from AM polarities.
Note I said "may"; my dummy example shows how the two can share the fundamental period, but this doesn't happen in OP's example. If there we set $f_m=100$, this is confirmed. Olli's answer fairly demonstrates how exact periodicity can be predicted from FT.


Why optical flow ? why not tracking I have a question that might be stupide !

Supposing that I try to detect object moving in a video or human action.

Many works are based on optical flow computation.

My question is why using OF is the best for motion detection. A tracking of pixels or region of interest will not be as informativa as OF .?

There are two versions of optical flow(OF): Feature based (sparse) or dense. In the dense version OF is applied to all the image pixels, while in the sparse one, only certain characteristic feature points are tracked. However, both approaches depend on the tracking of pixel quantities. This is fundamentally different than tracking the whole patch, because in return one obtains a full set of pixel-level correspondences. This is not just a vague statement, but rather carries the idea that the correspondence estimation is not constrained. Remember in the motion case, one would assume that the motion is "rigid" or "articulated" and imposes this prior into the tracking framework. If optical flow is used to estimate the dense trajectories, than this constraint is not assumed and one could as well track the deformable bodies, regardless of the deformation model.

Optical flow provides you more freedom and information about the tracked scene. Yet, my experience is that the motion tracking methods are more robust and reliable. But keep in mind that, by applying human models or temporal models on top of such tracked trajectories (from OF) researchers are capable of developing robust articulated / deformable / model-based tracking algorithms.

The math behind optical flow is the same, as in Lucas-Kanade tracker. The only difference, that tracking is usually applied to window or patch, and OF - for full image. So optical flow is motion tracking for whole image.


Guitar fundamental frequency vs pitch When a note is played on the guitar, is the fundamental frequency always the strongest one or could a harmonic be the strongest one (in terms of amplitude). When could that happen? Also, when playing a note of, for example, 440Hz, is 220Hz also a harmonic or only larger multiples?

When I produce a spectrogram of a B3 guitar note I get this:

As you can see, the frequency of B3 (246Hz) does not seem to be the strongest frequency. Its approximate half (near 130) seems to be the strongest one.

This is the audio file (the B3 note is clear).

To generate the spectrogram, I use a STFT taken from the Librosa library:

import plot as plt

def transcribe(filename):
  y, sr = librosa.load(filename, sr=40000)

  D = librosa.stft(y)
  plt.plot_spectrogram(D)

The plot.py module:

import matplotlib.pylab as plt

def plot_spectrogram(stft_output, xlabel='Time', ylabel='Hz'):
  plt.figure()
  librosa.display.specshow(librosa.amplitude_to_db(stft_output,
    ref=np.max), y_axis='log', x_axis='time')
  plt.title('Power Spectrogram')
  plt.colorbar(format='%+2.0f dB')
  plt.tight_layout()
  plt.savefig('static/plots/' + 'spectrogram.png')

Pitch is psychoacoustic.  Thus, a harmonic frequency can be stronger than the fundamental, without changing the pitch heard.  This is common in some stringed instruments.  And a note's pitch frequency could even be missing from a notes frequency spectrum.  For instance: An old telephone circuit can cut off frequencies below 200 Hz, but people still have no problem recognizing a low male voice pitched at around 100 Hz.

For many timbres, humans will hear the harmonic spacing as the note pitch.  In your plot, the harmonic spacing appears to be around 130 Hz, so that will be the note pitch perceived (if played at the same sample rate at the plot axis suggest).  Not 260.

The note/pitch is determined only by the first frequency peak on the DFT plot. All the other amplitudes are called harmonics that only add to the timber. These harmonics (timber) are what differences a guitar note, for example, from a note on a violin.


baffled by fft phase spectrum! A very simple MATLAB experiment:

f = 200;
fs = 1000;
t = 0: 1/fs : 1;
x = cos(2*pi*f*t);
plot(angle(fftshift(fft(x))));

And here's the output:

Now, made a minor change to the above code snippet; reducing the time duration by just 1 sample, as follows:

f = 200;
fs = 1000;
t = 0: 1/fs : 1 - 1/fs;
x = cos(2*pi*f*t);
plot(angle(fftshift(fft(x))));

And the phase spectrum goes totally crazy:

Questions:

In the first plot, I was hoping to see a zero phase at bin 700 which corresponds to the positive frequency of 200 in this example. That doesn't seem to be the case.
Secondly, I don't understand the linear parts of the graph in plot 1. I do appreciate phase components that could potentially exist due to the so-called numerical noise, but then how could that noise be so 'linear' in phase?
In the second plot, why removing just one sample would have such a drastic impact on the phase plot?
Am I doing something fundamentally wrong here?

You aren't doing anything wrong, but you also aren't thinking carefully about what you should expect to see, which is why you're surprised at the result. For question 1, your conjecture is close, but you actually have things backwards; it's numerical noise that's plaguing your second one, not your first one.

Pictures may help. Here's plots of the magnitude and phase for the first trial:

x = Cos[2.0 \[Pi] 200 Range[0, 1, 1/1000]];
fx = Fourier[x];
ListLinePlot[Abs[fx], PlotRange -> All]
ListLinePlot[Arg[fx], PlotRange -> All]

And the second one:

x = Cos[2.0 \[Pi] 200 Range[0, 1 - 1/1000, 1/1000]];
fx = Fourier[x];
ListLinePlot[Abs[fx], PlotRange -> All]
ListLinePlot[Arg[fx], PlotRange -> All]

So what's going on here? The second one is easiest to explain. First, the FFT of the second one has zero amplitude everywhere except at the two peaks visible in the magnitude spectrum; this follows since the definition of the FFT using 1,000 datapoints returns frequencies of the form $k/1000$ for $0\leq k\leq 999$, and hence your signal falls exactly on a frequency bin. As a result, at the other 998 points your signal is entirely due to machine noise due to floating point errors, and hence your phase spectrum is nonsense, because it's literally the phase of pseudorandom numbers.

For the first one, however, the definition of the FFT encompasses frequencies of the form $k/1001$ for $0\leq k \leq 1000$, whereas your signal frequency is $200/1000$, which is not of the form $k/1001$. As a result, your signal is broadened by spectral leakage, and will be nonzero almost everywhere. I will not comment on the physical form of the phase plot, but I will say that it admits a closed analytic form.

In general, I think graphs of phase angle alone are a really bad idea for conveying information, precisely for this reason; first, you can't tell if you're looking at the phase of low-amplitude garbage or actual signal, and second, it's not translation-invariant, and it's easy to get utterly bewildering graphs for simple inputs. Much better, if you're still looking for something which conveys phase information, is a graph which simultaneously portrays phase and amplitude information in the same visual manner, such as a plot where phase is encoded as hue and magnitude is encoded as brightness.

ADDENDUM: Here's a couple pictures from Mathematica which illustrate the principle I stated in the preceding paragraph:

hue = Compile[{{z, _Complex}}, {(1.0 Arg[-z] + \[Pi])/(2 \[Pi]),
Exp[1 - Max[Abs[z], 1]], Min[Abs[z], 1]},
CompilationTarget -> "C", RuntimeAttributes -> {Listable}];
L = 500;
data = Table[Boole[x <= 11 && y <= 11], {x, L}, {y, L}];
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {-5, -5}]], {L/2, L/2}]],
ColorSpace -> Hue, Magnification -> 1]
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {-4, -4}]], {L/2, L/2}]],
ColorSpace -> Hue, Magnification -> 1]
Image[hue@
RotateRight[
10 RotateRight[Fourier[RotateRight[data, {0, 0}]], {L/2, L/2}]],
ColorSpace -> Hue, Magnification -> 1]

All three images are 2D Fourier transforms of the same input signal (an $11\times 11$ square of 1's padded with zeros to a length of $500\times 500$), but the inputs have been cyclically rotated by 5, 4, and 0, and 200 datapoints. The magnitude spectra (encoded by pixel brightness) are identical, but the phase spectra are completely different! Phase encoding is done so that 1 maps to red, $i$ maps to green, $-1$ maps to cyan, and $-i$ maps to purple. This is what I mean when I say that phase spectra are non-shift invariant, and thus not amenable to human visual understanding. For example, with a cyclic shift of 200 datapoints, it's completely impossible to tell what's going on in the phase, since it just looks like static, but the input signal is no more complicated than the other input cases.

The Gaussian Waves site details the part about the phase and its random-like behaviour: just a matter of floating-point error as DumpsterDoofus said


Matlab version rpmordermap and rpmfreqmap Are those functions only available in Matlab 2015b version or they are included in 2015a version also? Thanks in advance.

At the very bottom of MATLAB manual pages there is always an info about the time when given function was first introduced. For both of rpm*map functions:

No, they are available from R2015b, as answered by @jojek.

Meanwhile, especially for older versions, you can find useful alternatives for Order-RPM function, for instance at MatlabCentral with Campbell diagrams, or Vold-Kalman order tracking code.


What is causing these spikes in the output of the low pass filter? I am low pass filtering the PCM of an arbitrary Bell 103 (300 bps) signal. However, the output of the low pass filter shows random spikes in amplitude.

I have tried to increase the bandwidth, but no luck.
Can someone explain why these spikes arise? My goal is to do some sort of envelope detection on this signal afterwards, to decode the 1's and 0's. For example by using a rectifier first.

The filter coefficients:

coeffs = [0.265842424474100 0.243276902464240 0.271096721997119 0.243276902464240 0.265842424474100]

The Matlab script used:

close all;
clear all;

% Generate a random bit sequence as input
bitsIn = randi([0 1], 16, 1)';

% Generate the amplitude of the FSK signal
amplitude = randi([7000 10000], 1, 1);

% Generate a randomized length of randomized startup noise
% x = randi([-20 20], 1, randi([0 3000], 1, 1));
x = []; % TODO

% FSK module the input bit stream
for i = 1:length(bitsIn)
    if (bitsIn(i) == 1)
        freq = 1270;
    else
        freq = 1070;
    end

    % Append the new PCM data for the new bit, with random phase
    phase = (2 * rand(1)) * pi;
    x = [x (amplitude * sin((2 * pi * freq * (0:26) / 8000) + phase))];
end

% Lowpass filter the resulting FSK samples
d = fdesign.lowpass('N,Fp,Fst,Ast',4,1070,1270,10,8000);
Hd = design(d,'equiripple');
y = filter(Hd,x);

figure();
plot(x);
title('Input data');
figure();
plot(y);
title('Low pass output');
```

I don't have MatLab with the DSP toolbox, but you're calling filter wrong.  Use:

y = filter(Hd, [1], x);

to set the feedback taps properly to $[1]$.

You should also inspect your filter design by using

freqz(Hd, [1], 512, 8000);

Although my hunch is that amplitude spikes in the output are directly caused by the filtering of the phase discontinuities you intentionally added when performing FM modulation.

Really, no modern equipment uses frequency switching to perform FSK.  FSK emitted from modern equipment is almost always continuous phase.  Are you really dealing with an ancient piece of Bell 103 modem hardware?

An example Octave script that generates the complex I/Q samples of continuous phase FSK is here:

https://dsp.stackexchange.com/a/53834/28112

if you need an example of how to generate CPFSK.


What exactly is the effect of a Hann Window on the FFT output? In the code I am trying to understand, there is an usage of Hann window multiple times (There's overlapping frames, the size of each is 512 and the overlapping factor is 75% so each time we capture  128 samples we build a new frame by shifting the previous one by 128 to the left). For each Frame we apply a Hann window before the FFT, so the first 128 samples were actually "framed" 4 times.

However I don't understand the purpose of this windowing. Can someone please tell me what its effect is on the FFT ?

See the following: http://www.dtic.mil/dtic/tr/fulltext/u2/a381695.pdf  especially p.631 Figs. 17+18.  As noted by the author, a 4 to 1 overlap using a Hann or Hamming window allows for proper reconstruction.

Also see: F. J. Harris, “On the use of Merged, Overlapped and Windowed FFT’s to Generate Synthetic Time Series Data with a Specified Power Spectrum,” Proc. IEEE 16th Asilomar Conf., 1982, pp. 316-321.

Convolution with the transform of a Von Hann window generates a less noisy looking FFT spectrum result than does convolution with a rectangular window.  One reason that this is so is because the Von Hann window attenuates the left and right sides of the data window so that there isn't a big difference between these two ends.  This makes the signal easier to represent with lower frequency basis vectors that are close to the same value at the two ends of the window, as are periodic-in-window sinusoids.


Histogram Equalization - mapping a certain portion of pixels to pure black and white How can I make sure to map, for example, 7% of the pixels to pure black (0) and pure white (255) in a grayscale image?

Compute the image histogram, then its running integral (i.e. an array T with as many entries as there are gray levels, such as T(k) = sum_{i=0, 1..k} (histogram(i)).

The last element of the integral is, of course, the area of the image.
For the lower threshold, start from the left and find the level at which T exceeds your desired fraction. Analogously, for the high threshold, start from the right and find the level at which T goes below your high fraction.

Although your claim does not make much sense to me, you may do following steps to ensure get what you described in your question.

Construct some discriminative feature space based on the original image.
Find two pixels that are separated most in this feature space.
Pick one pixel as your landmark, and sort all pixels with respect to their distance away from this landmark.
Set the first PERC% of these pixels to black, and the reset to white.

The "discriminative feature" can be any feature. The whole point here is to make different pixels different values. You may directly use raw pixel intensity if you want, but keep in mind of possible cases that there are 10% of 0 pixels, but you only want to set 7% of all pixels to black.

The general idea here is to convert all pixels in an image into a new "ranked" feature space. In this way, selecting PERC% of pixels is equivalent to set a threshold of value PERC to this feature map.


How do I analyze a music file ( only contains 5 seconds of note A) using the FFT plot to determine the note being played? To perform this task why can't I just import the file, and directly perform the fft and use the stem plot to determine the note? From the index k of the peak I can do w = 2pi*k/N where N is total number of points.

Fs, x = wavfile.read('piano.wav')
X = np.fft.fft(x)
plt.plot(abs(X))
plt.show()

However I'm not sure about the correct way to map the fourier coefficients from the plot to the frequency to determine the note.

Please advise.
Thanks.

An FFT magnitude peak frequency (derived from a single FFT result index) is not reliably the same thing as the music pitch frequency or the musical note.  Note pitch is a psychoacoustic phenomena, and more reliably related to periodicity than the (sinusoidal) spectral peak frequency (e.g. a periodic/repeating waveform shape does not need to look like a sine wave).  See FFT Pitch Detection methods: Autocorrelation or other? for a list of some better methods for estimating note pitch frequency or the musical note.

this was taken from this article here: http://pythonforengineers.com/audio-and-digital-signal-processingdsp-in-python/

frame_rate = 48000.0
infile = "test.wav"
num_samples = 48000
wav_file = wave.open(infile, 'r')
data = wav_file.readframes(num_samples)
wav_file.close()
data = struct.unpack('{n}h'.format(n=num_samples), data)
data = np.array(data)
data_fft = np.fft.fft(data)
# This will give us the frequency we want
frequencies = np.abs(data_fft)
print("The frequency is {} Hz".format(np.argmax(frequencies)))


Real-time frequency recognition in python I'm a newcomer to signal processing, so forgive my ignorance. I'm working on a project in Python that will play a .wav file and pick out the dominant frequency as it plays. I've found a few examples on the web, but none seem to be accurate for stereo files. It seems like I need to do something with an FFT, but I don't know too much about it.
Any help would be greatly appreciated. Thanks!

Here's a recent blog post that may be helpful:

https://jakevdp.github.io/blog/2013/08/28/understanding-the-fft/

You want to install numpy (and probably scipy, too), and read up on the Discrete Fourier Transform (DFT; this is what the FFT computes).  I'm not that up on Python media packages for the playback part. Python's wave module may be enough.

The algorithm could be as simple as computing the FFT of the last, say, 4096 samples, getting the index of the maximum absolute value, and multiplying by the bin width.  Or you can get fancier and use a Hanning window, peak interpolation, overlap processing, and exponential averaging.


Can I see from the PSD if a system is a LTI system? I am studying the meaning of the power spectral density (PSD) and I am wondering whether I can see from the PSD of the input signal and output signal if a system is a LTI system?

Suppose that the input signal to the system has power spectral density (PSD)
$S_i(f)$ while the output has PSD $S_o(f)$. Then, if the system is an LTI
system, then it must be that
$$S_o(f) = |H(f)|^2 S_i(f) \tag{1}$$
where $H(f)$ is the transfer function of the system. So, for a given pair
of PSDs, is there a $H(f)$ such that $(1)$ holds? One obvious case when
$(1)$ cannot hold for any choice of $H(f)$
is when the supports of $S_i(f)$ and $S_o(f)$ are different (in the sense
that there are frequencies $f$ for which $S_o(f) > 0$ while $S_i(f) = 0$).
So in this
case we can say unequivocally that the system is not an LTI system.

What about when the supports of $S_i(f)$ and $S_o(f)$ are identical or
the support of $S_o(f)$ is a subset of the support of $S_i(f)$ (think
bandstop or notch filter)?
Well, in this case, we do have a $|H(f)|^2$ from which we can
deduce $|H(f)|$. Throw in a conjugate-symmetric $\angle H(f)$ and we
have constructed the transfer function $H(f) = |H(f)|e^{\angle H(f)}$
of an LTI system whose output
PSD is $S_o(f)$ when the input PSD is $S_i(f)$. Note that there can
be infinitely many LTI systems that will provide this specific
input-output relationship.

But, the above does not allow us to say that the system that actually
produced the given output PSD from the given input PSD is an LTI
system, let alone one having the specific transfer function that we have
constructed. The actual system that we are looking at might be nonlinear
or time-varying, etc. We cannot say. After all, the
$|H(f)|e^{\angle H(f)}$ that we have come up with
might be a physically unrealizable transfer
function.

Yes, the PSD of a time signal can give you infos about the linearity of the system under study.

For example, assume you are studying a symmetric duffing's oscillator described by the equation:

$m \ddot{x} + c \dot{x} + k x + k_2 x^3 = p(t) $

$p(t) = A \sin{\omega t}$

If you compute the PSD of the acceleration signal you will see a peak at the excitation frequency and then some harmonics generated by the non linearity.

Since the nonlinearity is odd, the harmonics will be at odd multiples of the excitation frequency, so if $\omega$ is the excitation frequency, the harmonics will be at $3\omega$,  $5\omega$ and so forth.

For linear system these harmonics are not present, because the response of a linear system excited by a monoharmonic signal will be monoharmonic at the same frequency (after transients have died out).

I suggest you to read the book written by Worden and Tomlinson if you are interested in nonlinear dynamics.


LPF - signal values unaffected at specific times Is it possible to design an LPF that has an output identical to the input at specific points in time domain (the rest of the input waveform can get filtered/distorted)?
Is there a general name/technique for this kind of thing (assuming it is possible), so that I can search for more information on this topic?

What you probably are looking for are Nyquist/$M$th-band filters. The impulse response of these have the following property (assuming a non-causal impulse response centered around tap 0 for ease of exposition)

$h_n = \left\{ \begin{matrix} \frac{1}{M} & n=0\\ 0 & n = kM, k=\pm 1, \pm 2, \dots\end{matrix} \right.$

This means that every $M$th sample is untouched (the $1/M$ term is for scaling purposes). Often it is said that these filters have zero inter symbol interference, assuming that a symbol arrived every $M$th sample.


Intuition behind the Continuous Wavelet Transform? I was thinking sometime back about how to explain the Continuous Wavelet Transform ELI5. So this is what I came across.

The correlation of two exact signals is 1. So if I have an input signal $f(x)$ made up of an array of frequencies, how can I find out what frequencies exist at what points? Well, slide a signal $m(y)$ where $-\infty > y > +\infty$ over $f(x)$, and at those points where the co-relation of these signals is 1, well, those frequencies are present, at those times. This if of course a Continuous Wavelet Transform. Am I correct?

I think that the best way to explain the CWT is to start by explaining the Fourier Transform, then move on to explaining the Short-Time Fourier Transform, and then finally explain the CWT as a variation of the STFT.

The Fourier Transform exploits the fact that any decently behaved function can be represented as a sum of sinusoids (i.e. a Fourier series) and that the sinusoid basis posses the property of orthogonality:

$$ \int_{-\infty}^{\infty}  \sin(nx)\,\sin(mx)\,dx = \int_{-\infty}^{\infty}  \cos(nx)\,\cos(mx)\,dx = \begin{cases}
1,  & \text{if $n=m$} \\
0, & \text{if $n\neq m$}
\end{cases}$$

So, since:

$$ e^{iax} = \cos(ax)+i\,\sin(ax)$$

the Fourier Transform is simply doing this integration for all frequencies and keeping track of which outputs are zero (i.e. that frequency is not in the signal) and which are non-zeros (i.e. that frequency is in the data and its output is scaled by how much of it is in there):

$$ S(f) = \int_{-\infty}^{\infty} s(t)\,e^{i2\pi ft}dt $$

In this case you are doing this integration over the entire signal so you can't really tell if the frequency content is changing from the beginning of the signal to the end. One way around this is to compute the Short-Time Fourier Transform: i.e. window the signal, calculate the Fourier transform of the windowed signal, store it, then shift the window down a bit and repeat for all shifts:

$$ S(\tau,f) = \int_{-\infty}^{\infty} w(t- \tau)\,s(t)\,e^{i2\pi ft}dt  $$

where

$$ w(\tau,t)=\begin{cases}
1,  & \text{if $\tau\approx t$} \\
0, & \text{if $\tau \not\approx t$}
\end{cases} $$

The key thing here is that you are calculating the typical Fourier transform but of a new signal that only exists in a localized part of the t-axis. To emphasize this, you can see the new signal whose Fourier transform we are calculating by associating:

$$ S(\tau,f) = \int_{-\infty}^{\infty} [w(t-\tau)\,s(t)]\,e^{i2\pi ft}dt  $$

And here's a graphical example of this showing the new signal for different $\tau$ values and on the right are a few sinusoids to represent what we are using to decompose the signals (i.e. our basis, or kernel).

But, we can also change the association as such without changing the outcome:

$$ S(\tau,f) = \int_{-\infty}^{\infty} s(t)\,[w(t-\tau)\,e^{i2\pi ft}]dt  $$

So this means that instead of windowing our signal, we are windowing our basis functions. But here's the kicker, if we are windowing our basis functions, we don't have to use a constant-size window since we know that a basis function of high frequency will need a shorter window than a basis function of low frequency. This is the whole point of the CWT. It is a decomposition of a signal by "wavelets" (i.e. windowed sinusoids in this case) where the windowing is adaptive to the sinusoid frequency. If the we choose a Gaussian window (as I have chosen in these examples), then our wavelets are called Morlet wavelets (or Gabor wavelets, in some literature).

Finally, you can generalize this for any choice of wavelet that you want. In that generalization, you can describe your wavelet basis functions as being "stretched" and "squeezed" versions of some arbitrary "mother" wavelet, $ \psi $. And so the previous equation can now be written as the final form of the CWT:

$$ S(a,b) = \frac{1}{\sqrt b}\int_{-\infty}^{\infty} s(t)\,\psi (\frac{t-a}{b})dt $$

where $a$ is now what used to be our $\tau$ (i.e. a time shift), and $b$ is called the "scale" which is just a parameter to stretch and squeeze the wavelet (similar to what our parameter $f$ was except now the interpretation is more difficult). And the only reason you have the $\frac{1}{\sqrt b}$ up front is to normalize the wavelets so that they have the same "energy" and you end up comparing apples to apples in your time-frequency representation.

I hope this helps!

-Antonio

You might be interested in the following link where wavelet transform is intuitively presented. It starts from Fourier transform and points out the problem of loss of localisation in frequency domain. It is followed by Short Time Fourier Transform (STFT). Duality in time frequency is also stated for STFT and then the concept of SCALE is introduced followed by THE WAVELET TRANSFORM

https://www.youtube.com/watch?v=qy6r5TE7YxU&index=10&list=PL_1-_Q4VI1Sdfh6hLIlDnZj6TY0pXQwkN


How do I optimize the window lengths in STFT? I have many EEG signals and I want to analyze them using linear methods such as STFT (Short Time Fourier Transform). In STFT , How can I optimize the analysis window length, to reflect the frequency spectrum of each analysis window in a proper way?

This is the classic "uncertainty principle" of the Fourier Transform. You can either have high resolution in time or high resolution in frequency but not both at the same time. The window lengths allow you to trade off between the two.

If you want to detect "events" in your EEG signal with a resolution of say 10ms, then this should be your window length. This will give you a frequency resolution of about 100 Hz.

Usually wide window size gives better frequency resolution but poor time resolution and vice-versa. Look at this example where I generated a spectrogram of a sine wave with 5kHz and sample rate of 22050Hz, from my C++ code.

The above spectrogram has window size of 2048 samples and overlap of 1024 samples.

Look at this spectrogram:

This one has window size of 512 samples and overlap of 256 samples.

Can you see the difference? The first one has better frequency resolution than the second one. But the second one has better time resolution when compared to the first one. So, choosing window size depends upon your application. If you are dealing with speech samples to track pitch, choosing larger window size should be the appropriate one.


MMSE- Channel estimated based - DFE Equalization For HF ionospheric channels ( frequency range 3 up to 30 MHz), I Want to make channel equalization in time domain, using MMSE criterion. I have to do it using the channel estimation.

I have some questions:

Can some one elaborate more on the closed-form solution solution for the feedforward and feedback FIR filters ?
And what is the exact relation between the number of taps of the filters and the channel delay spread ?

Any help will be really appreciated

Thank you

It seems that this PhD dissertation covers part of your question:
« Improved receivers for digital High Frequency communications: Iterative channel estimation, equalization, and decoding (adaptive turbo equalization)»
https://ntnuopen.ntnu.no/ntnu-xmlui/bitstream/handle/11250/249838/126139_FULLTEXT01.pdf?sequence=1&isAllowed=y

Edit:
4.2 tapped delay line model, 4.11-4.13 relates the (complex baseband) tapped delay line to physical channel multipath. (p.p 57)


3D reconstruction of an image based on depth of intensity levels? Is there anyway of segmenting objects based on how far they are from the viewer?

Can color values assess this sort of a thing for us?
How will intensity levels help in determining how far the object is from the viewer?

Another image:

This is an extremely difficult task, one which is a very active line of research. I've managed to find a semi-recent paper on the subject. I won't go into the details, but here's a few things that this paper found that can be used.

Textures: If you see a given texture, they can help you figure out how far something is away. The wood texture here would be quite ideal for this.
Fixed object sizes: Identify common objects, see how large they are, and use them to help with depth analysis. The phone book, soda bottle, etc should work well with this.
Focus: If you know how the picture was take, focus can help you identify ranges.
Shadows- If you know where the light was coming from, shadows can help you identify how tall objects are, and thus how far they are away.

See this code from matlab website: "Depth from gradient":

http://www.mathworks.com/matlabcentral/fileexchange/45269-dfgbox


How to do the highest resolution FFT? I wish to do the FFT of streaming audio data coming into a computer. I am just interested in frequencies from 80 Hz up to 1000 Hz but with a frequency resolution of 3 to 4 Hz and no more.

Basically my question is , as my sampling frequency (Fs) is 44100, will my FFT window size have to be no less than 11025 (N) ?  (because spacing of FFT bins = Fs/N). Does this mean that I am completely limited to taking windows of data any higher than 4 times a second with my desired resolution.

Unfortunately I need to be able to capture windows of data 8 times a second and also have a resolution of 3 to 4 Hz.

Is this actually possible? Many thanks in advance.

The FFT size can be longer than the window size.  For instance you can use audio sample windows of 5512 (about 8 per second) and zero-pad them to length 16384 for each FFT.  If noise level is low enough and the separation between frequency peaks is great enough, the interpolation provided by this zero-padding can provide better frequency estimation and plot resolution.

The "highest resolution" FFT reasonable depends on the signal-to-noise ratio as well as the window length, as after a certain point, interpolation and frequency estimation methods will be affected as much by nearby spectral noise as by the spectral peaks that one is trying to estimate at "high resolution".

You can zero pad your "windows" of data to increase the frequency resolution. to get 3-4 Hz resolution, you would need to take a 16384 point FFT (gets you about 2.5Hz resolution). you could then take 1/8th of your sampling rate in samples (5513 samples) and add zeros to the end until you get 16384 samples. By taking the FFT of this, new block with zeros at the end, you still get a 16384 point FFT, which gives you the frequency resolution you need, but still be able do your fft 8 times a second.

Please note that you need a number of samples that is a power of 2 to do an FFT, If you are doing the fft in matlab, i believe it will add zeros to the end of your signal if you don't give it enough samples. If you can get by with a resolution of 5.4Hz, then you can do a 8192 point FFT.


How do I set window length and shift length for computing STFT on Scipy? Does nperseg=32, noverlap=16 mean the window length is 32mS with a shift of 16mS?
Refer here for more details

for this function:

nperseg will set the number of samples to output from your STFT.
noverlap will specify how many old samples you want to take to calculate your future STFT.

Let's say you want a stft of size 256 with an overlap of 50% you will set noverlap = 128.
128 samples will then be taken from the old segment that was used to calculate the previous STFT...
For more information about STFT and the interest of overlap I invite you to watch this discussion.
STFT: why overlapping the window?


Does the instantaneous frequency contradict the uncertainty principle? The uncertainty principle states that there is a trade off between time and frequency. So, finding frequency components at specific time is impossible. However, the instantaneous frequency measure the frequency as a function of time. Which means using the instantaneous frequency, the frequency components could be found for a signal at a specific time. How can you interpret this? Why don't we use the instantaneous frequency for time-frequency analysis?

The uncertainty principle works in the presence of (an uncertain amount of) noise or other signals (including possible harmonics), corrupting the exact phase, and thus the rate of change of phase of the signal of interest.  If the phase is corrupt, or mixed with the phase of other signals, then deriving an instantaneous frequency from the 1st derivative of that phase might produce nonsense.  Time-frequency analysis might be one way to (statistically?) separate information about the signal of interest out of these potentially existing "corrupting" influences.

Whereas the phase of a perfectly analytic signal including zero additive noise is better defined.


Is it a good idea to whiten a colored noise in order to apply formulas for AWGN? I have posted a question before regarding the detection of a signal under colored Gaussian noise:
Energy detection in presence of colored Gaussian noise

Theoretical formulas for setting the detection threshold are widely used for AWGN.

Is it a good idea to just perform a pre-whitening to the colored noise and then apply the well-know formulas for energy detection under white Gaussian noise?

Do you have any happy experience using this approach in practice?

The short answer is yes, whitening noise is usually recomended (if you know the spectrum of the noise) but whitening also whitens your signal.  Detecting random signals in colored noise is a well known problem is sonar signal processing.  There are a number of derivation that I like particularly:

Nielsen, Richard O. Sonar signal processing. Artech House, Inc., 1991.

and Burdic

@book{RN17,
   author = {Burdic, William S.},
   title = {Underwater acoustic system analysis},
   publisher = {Prentice Hall},
   ISBN = {0139476075},
   year = {1991},
   type = {Book}
}

Also the updated

@book{van2004detection,
  title={Detection, estimation, and modulation theory, part I: detection, estimation, and linear modulation theory},
  author={Van Trees, Harry L},
  year={2004},
  publisher={John Wiley \& Sons}
}

in section 4.3

and if you can find it, in the old volume 3

If you look, a copy of Eckart's  classic paper can be found at Scripts.

Assuming we are talking about a linear system (i.e. within its dynamic range its behavior does not change with the strength of the incoming energy), then it is never worth adding extra noise energy in hopes of improving signal detection. Consider that after removing any out-of-band noise with a simple filter, any reasonable technique beyond that must operate on in-band noise. A short snapshot of in-band noise is almost always colored, anyway. So a "good" AWGN technique must, as a rule, also be a good colored noise technique.


Sample and reconstruct a real exponential (just one period) I have a function with an equation:

$$C = 1.6925\left( e^{-0.136t}-e^{-1.192t}\right) $$

Where $C$ is real and $t$ represents time in hours. Beneath is the representation of my function.

I am trying to find the best sampling step (based on Nyquist-Shannon theorem) to sample and reconstruct this function using a Fourier Transform.

After using a Fast Fourier Transform algorithm in MATLAB (fft), by setting:

x= fft(c);
xmag=abs(x); (amplitude estimation)

I find the below spectrum for the function for a $t$ between $0$ and $100$ sampled at $1000 \text{Hz}$  , with the highest amplitude at:

bin(1) = 1.1025e+04

and the second highest at bin(100001) with an amplitude of 9.9946e+03. the bin number 100001 represents also the length of xmag.

So my question is how would you interpret the spectrum? Specifically I would like to know how to determine the highest frequency at which I can sample this function and reconstruct it without aliasing?

First, since $t>0$ in your case, you can write your function as

\begin{equation}
C(t)=1.6925\left(\exp^{-0.136t}- \exp^{-1.192t} \right) u(t)
\end{equation}
where $u(t)$ is a unit-step function. Then, denote the FT of $C(t)$ as $C(F)$, which is given by
\begin{equation}
C(f)=1.6925\left(\frac{1}{0.136+2\pi jf}- \frac{1}{1.192+ 2\pi jf} \right)
\end{equation}

The spectrum $C(f)$ is plotted below. We can see that $C(f)$ is almost zero after 5 Hz. You can choose the sampling frequency according to the Nyquist criterion as $f_{\max} \geq 10$ Hz. You have taken a very high sampling rate (1000 Hz) that is 100 times more than the required rate. When I took 10 Hz to sample the signal above and took the FFT I got the results as follows  (in order time domain $C(t)$, frequency domain signal $C(f)$, frequency domain domain signal obtained via FFT, that is, FFT(C(t)). Note the inherent periodicity of the FFT spectrum)

And here is the MATLAB/Octave code for your reference..

clear all; clc;close all;

%% Time-domain signal
Ts=0.1; % Sampling time (obtained from the frequency domain plot)
t=0:Ts:60; % in sec
C=1.6925* (exp(-0.136*t)-exp(-1.192*t));
figure,plot(t,C);
title('Time domain signal');
xlabel('Time in sec')

%Frequency domain signal
f=0:.1:10;
F1=1.6925*((1./(0.136+(j*2*pi*f)))-(1./(1.192+(j*2*pi*f))));
figure,plot(f,abs(F1));
title('Frequency domain signal');
xlabel('Frequency in Hz')
% Via FFT
F2=fft(C);
L=length(F2);
figure,plot((0:L-1)*(1/(L*Ts)),abs(F2));
title('Frequency domain signal via FFT of C(t)');
xlabel('Frequency in Hz')

How did you arrive at sampling frequency of 1000 Hz?  Ideally speaking, you should find the Fourier transform of the given function analytically.  This can be done using the FT tables for exponentials.  Then you should find the highest frequency beyond which the FT is almost zero i.e less than the smallest number represented by finite precision you are using. Twice or more of this frequency is your sampling rate.  How many samples should you take?  The signal is almost zero beyond 5 x 10$^4$  as can be seen from the figure.  So you take samples upto this time point.  Now you take the FFT of these samples and it would be equal to the FT you analytically calculated i.e it would be sampled version of the FT.  If you used sampling rate less than this procedure explained, the resulting FFT would be aliased version.


Filtering out a specific frequency in an analog signal I'm wondering how I would go about filtering out a regular sine-wave signal with a constant frequency and voltage.

I'm really new to the whole DSP-area but imagine I would have a regular sinewave at a specific frequency (let's say 30Hz), and on top of that signal I have a bunch of other small peaks (which are the ones I would want to sample and read). What method would I choose?

Here is a visualization of what I meant:

Excuse my terrible Paint-skills but how would I filter out that sinewave, and just keep the small spikes in voltage? What might be worth sharing is that the time between the spikes does have an impact on my measurement of the spikes, or to rephrase: The amount of spikes within x-units of time does have an impact.

Any guidance to the right direction would really be appreciated.

Thank you.

A notch filter is a good approach in general for removing an unwanted frequency. But if you are certain of the exact frequency, phase and amplitude of the unwanted sine wave, then you may also consider adding the noisy signal to a sine wave of the same frequency and amplitude but 180 degrees out of phase with the unwanted sine wave.

To illustrate the notch filter example, the plot below shows an attempt to generate a sine wave with some odd impulses (black), and then filter it with a notch filter knowing the exact frequency (red).

R code for implementing it below.

# 24782
# install.packages('signal')

T <- 1024
t <- 0:(T-1)

fs <- 1000
omega <- 2*pi*30/1000
phi <- 0.9279835

x <- sin(omega*t + phi)

changed <- sample(1:T, 10)

y <- x
y[changed] <- y[changed] + 3

alpha <- 0.9
num <- c(1, -2*cos(omega), 1)
den <- c(1, -2*alpha*cos(omega), alpha*alpha)

yf <- signal::filter(num, den, y)

plot(t,y, type="l", lwd=10)
lines(t,yf, col="red", lwd=4)

legend(850, 4, c("Original", "Filtered"), lwd=c(2.5,2.5),col=c("black","red"))

You can sample the raw signal and then use a bandstop filter to filter out the 30 Hz frequency component.

http://www.mathworks.com/help/signal/examples/introduction-to-the-filter-design-and-analysis-tool-fdatool.html

Fdatool is a MATLAB tool simply to use


GPS CA Signal Acquisition The coarse acquisition C/A GPS signal has a symbol rate or chip rate of 1.023 Mcps with code period of 1 ms. The code bandwidth is +/-1.023MHz (to 1st sinc nulls). I'd like to record this signal using a USRP.

What's the minimum complex sample rate that I'd need to use so I could acquire this signal in software? 1 sample/chip = 1.023e6? 2 samples/chip = 2.046e6?

Since the signal is coming from a satellite, there will be an unknown carrier frequency offset due to Doppler. From what I've read, for GPS the Doppler is +/- 10 KHz. What's the maximum residual frequency offset allowable when doing the despreading of the CA code? 500 Hz? Or, what's the SNR loss due as a function of frequency missmatch due to residual Doppler?

What's the maximum residual frequency offset allowable when doing the despreading of the CA code? 500 Hz? Or, what's the SNR loss due as a function of frequency missmatch due to residual Doppler?

The maximum frequency offset is dependent on SNR required for acquisition as the roll-off of correlation versus frequency offset is due entirely to the duration of the correlation.

The magnitude of correlation vs frequency offset is a Sinc function with the first null at 1/T Hz where T is the duration of the correlation interval in seconds.  So if your correlate over 1 full PRN period, which is 1 ms, the correlation will go to zero at a 1 KHz Doppler offset.  Note the the PRN is repeated 20 times before a data transition (the data rate if 50 Hz), so you could actually correlate over 20 ms for 13 dB more processing gain (Processing gain is 10Log(N)), at the expense of reducing the Doppler tolerance (the first null would appear with a 50 Hz offset in this case!).

A good demonstration of this is the DFT itself for those familiar with it. Each bin in the DFT is a correlation over the time interval of the sequence to the particular frequency in that bin. The DFT can be described as a bank of filters with each filter approaching a Sinc magnitude response as N gets larger (when a rectangular window is used), with the first null appearing in the adjacent bin and then additional nulls in every other bin (hence we show that each bin is uncorrelated to the others). Given the time domain sequence is N samples long in time, and there are N bins also in frequency, then the spacing between bins is 1/N, as we would expect by the description above.

This is demonstrated in the generic graphic below showing the correlation to a 3 Hz sine wave (red waveform) buried in noise over a 2 second correlation interval.

From this we also see how the post correlation magnitude is reduced as a function of frequency offset, while the noise, if white, would not be affected- from which we can determine SNR vs frequency offset.

We also see the trade space involved between processing gain and Doppler offset tolerance. If the received SNR is strong, we can reduce the duration for the correlation and then benefit from a wider frequency acquisition window. For example with GPS we could correlate over half a PRN sequence (0.5 ms) at the expense of 3 dB in processing gain, with the benefit of having the first null out at 2 KHz instead of 1 KHz.  In practice I have typically used the metric of half the main lobe as my Doppler range for the correlation (so for GPS would be 1 KHz wide or +/- 500 Hz when correlating over 1 ms), but in weaker signal conditions a smaller search range may be necessary in order to maximize the SNR required for acquisition.   (Meaning how far I would step the frequency to repeat the correlation again while searching for signals; typically with a 1ms correlation, I would step 1 KHz but if I felt more SNR may be necessary I may search in finer steps, I could also increase the correlation time which would also require finer frequency steps in the search but here I was just referring to being less sensitive to the roll-off within a set correlation time by not stepping as far in frequency before repeating the correlation again).

Also see this post related to GPS acquisition and a joint delay and frequency acquisition approach that can acquire in one PRN sequence (at the expense of significant processing): GPS signal acquisition

Also to note, once within a Doppler bin, the I and Q correlated outputs can then used to determine the precise Doppler offset and do carrier tracking with implementations described here (where we see the algorithm to determine the precise frequency offset from two adjacent I Q correlation outputs is quite simple):
High modulation index PSK - carrier recovery

The critical sampling rate is the bandwidth with complex baseband sampling, where a single sample is in fact a complex number.

So, in theory, 1.023 MS/s would be enough. You'd typically go higher to mitigate the filtering imperfection effects as well, as you've noticed, problems due to doppler. (

Also, not all devices can sample at all rates, and 1.023 MHz isn't really an easy rate to achieve.

Look at existing SDR GPS receivers like GNSS-SDR to see how they sample.

  What's the maximum residual frequency offset allowable when doing the despreading of the CA code?

That's a bit of a different question.

Let's do a quick estimate here:

If I remember correctly, the spreading sequence is 10230 chips in length
That means a symbol takes 10230/1.023e6 = 10⁻² s
Modulation is BPSK -> maximum allowable phase shift $\frac \pi2$, i.e. ¼ of a cycle
$\frac{\frac 14}{10^{-2}\,\text{s}}= 25\,\text{Hz}$ is the maximum allowable residual frequency error.


Calculating the inverse filter for the (exponential) sine sweep Method I am currently working on my Bachelor-Thesis about Real-Time convolution and Impulse Response Measurements. After reading different papers about the (exponential)sine sweep method I didn't find an answer on how to calculate the inverse filter for deconvolving the impulse response.

In the papers I read, it is described as the Time Reversal Mirror and after that some kind of scaling.

Could someone explain, how to calculate the inverse filter for a given sine sweep?
If you could add an example, or an algorithm it would be appreciated

Assuming that your Exponential Sweep Sine was generated using the formula:

$$x(t)=\sin\left(\frac{2\pi f_1 T}{R}\left(e^{\frac{t R}{T}} -1\right) \right)$$

where:

$f_1, f_2$ - Initial and final frequency of the sweep

$T$ - Duration of the sweep

$R = \ln\left(\frac{f_2}{f_1} \right)$ - Sweep rate

Then the inverse filter is calculated by scaling the amplitude of time reversed $x(t)$ by:

$k = e^{\frac{tR}{T}}$

Which will result in an exponentially decaying sweep:

$f(t) = \frac{x_{inv}(t)}{k}$

Example in Python:

#!/usr/bin/env python

from __future__ import division
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

def dbfft(x, fs, win=None):
    N = len(x)  # Length of input sequence

    if win is None:
        win = np.ones(x.shape)
    if len(x) != len(win):
            raise ValueError('Signal and window must be of the same length')
    x = x * win

    # Calculate real FFT and frequency vector
    sp = np.fft.rfft(x)
    freq = np.arange((N / 2) + 1) / (float(N) / fs)

    # Scale the magnitude of FFT by window and factor of 2,
    # because we are using half of FFT spectrum.
    s_mag = np.abs(sp) * 2 / np.sum(win)

    # Convert to dBFS
    ref = s_mag.max()
    s_dbfs = 20 * np.log10(s_mag/ref)

    return freq, s_dbfs

if __name__ == "__main__":
    # Sweep Parameters
    f1 = 10
    f2 = 100
    T = 3
    fs = 1000
    t = np.arange(0, T*fs)/fs
    R = np.log(f2/f1)

    # ESS generation
    x = np.sin((2*np.pi*f1*T/R)*(np.exp(t*R/T)-1))
    # Inverse filter
    k = np.exp(t*R/T)
    f = x[::-1]/k
    # Impulse response
    ir = sig.fftconvolve(x, f, mode='same')

    # Get spectra of all signals
    freq, Xdb = dbfft(x, fs)
    freq, Fdb = dbfft(f, fs)
    freq, IRdb = dbfft(ir, fs)

    plt.figure()
    plt.subplot(3,1,1)
    plt.grid()
    plt.plot(t, x)
    plt.title('ESS')
    plt.subplot(3,1,2)
    plt.grid()
    plt.plot(t, f)
    plt.title('Inverse filter')
    plt.subplot(3,1,3)
    plt.grid()
    plt.plot(t, ir)
    plt.title('Impulse response')

    plt.figure()
    plt.grid()
    plt.semilogx(freq, Xdb, label='ESS')
    plt.semilogx(freq, Fdb, label='Inverse filter')
    plt.semilogx(freq, IRdb, label='IR')
    plt.title('Spectrum')
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Amplitude [dBFS]')
    plt.legend()

    plt.show()

And output:

For Robert, here is the spectrum plot for case without amplitude modulation of the inverse filter:

Relevant literature:

  Q. Meng - Impulse Response Measurement with Sine Sweeps and Amplitude
  Modulation
  Schemes

  A. Novak - Nonlinear System Identification Using Exponential
  Swept-Sine
  Signal

  K. Vetter - ExpoChirpToolbox - a Pure Data implementation of ESS
  impulse response
  measurement

@jojek answer is very good.

Long history short, the inverse ESS should be the ESS time reversed with amplitude scaling so after convolving with the ESS the frequency is constant in the [f1 f2] range.
In other words the closest to a diraq impulse.

The matlab code would be (approximately):

invSweep(n) = sweep(N-1-n) .* (f2/f1)^(-n/(N-1))

invSweep(n) = flip(sweep) .* (f2/f1)^(-n/(N-1))

N being the duration in samples of the ESS (sweep)

Perhaps there's some error in the style N instead of N-1.


Radix-4 FFT versus Radix-2 Is a radix-4 implementation faster than a equivalently well coded radix-2 FFT?  And if so, why would it be faster?

It depends. Theoretically you can save a few multiplies with a radix-4 as radix-4 has a 1/4th the number of butterflies and 3 mpy + 8 adds per butterfly (if properly structured) and the radix 2 has 1 mpy + 2 adds per butterfly.

So in terms of multiplies it's a bit better, however there is higher complexity in terms of code structure, exception handling, coefficient management, register management, digit-reverse addressing, etc.

So it's only an advantage if the number of mpy is the limiting factor which for most hardware these days is not the case.

In radix 2, the number of sample is in terms of power of 2 power but in radix 4 the number of samples belong is a power of 4.