Why so many methods of computing PSD? Welch's method has been my go-to algorithm for computing power spectral density (PSD) of evenly-sampled timeseries.  I noticed that there are many other methods for computing PSD.  For example, in Matlab I see:

PSD using Burg method
PSD using covariance method
PSD using periodogram
PSD using modified covariance method
PSD using multitaper method (MTM)
PSD using Welch's method
PSD using Yule-Walker AR method
Spectrogram using short-time Fourier transform
Spectral estimation

What are the advantages of these various methods?  As a practical question, when would I want to use something other than Welch's method?

I have no familiarity with the Multitaper method. That said, you've asked quite a question. In pursuit of my MSEE degree, I took an entire course that covered PSD estimation. The course covered all of what you listed (with exception to the Multitaper method), and also subspace methods. Even this only covers some of the main ideas, and there are many methods stemming from these concepts.

For starters, there are two main methods of power spectral density estimation: non-parametric and parametric.

Non-parametric methods are used when little is known about the signal ahead of time. They typically have less computational complexity than parametric models. Methods in this group are further divided into two categories: periodograms and correlograms. Periodograms are also sometimes referred to as direct methods, as they result in a direct transformation of the data . These include the sample spectrum, Bartlett's method, Welch's method, and the Daniell Periodogram. Correlograms are sometimes referred to as indirect methods, as they exploit the Wiener-Khinchin theorem. Therefore these methods are based on taking the Fourier transform of some sort of estimate of the autocorrelation sequence. Because of the high amount of variance associated with higher order lags (due to a small amount of data samples used in the correlations), windowing is used. The Blackman-Tukey method generalizes the correlogram methods.

Parametric methods typically assume some sort of signal model prior to calculation of the power spectral density estimate. Therefore, it is assumed that some knowledge of the signal is known ahead of time. There are two main parametric method categories: autoregressive methods and subspace methods.

Autoregressive methods assume that the signal can be modeled as the output of an autoregressive filter (such as an IIR filter) driven by a white noise sequence. Therefore all of these methods attempt to solve for the IIR coefficients, whereby the resulting power spectral density is easily calculated. The model order (or number of taps), however, must be determined. If the model order is too small, the spectrum will be highly smoothed, and lack resolution. If the model order is too high, false peaks from an abundant amount of poles begin to appear. If the signal may be modeled by an AR process of model 'p', then the output of the filter of order >= p driven by the signal will produce white noise. There are hundreds of metrics for model order selection. Note that these methods are excellent for high-to-moderate SNR, narrowband signals. The former is because the model breaks down in significant noise, and is better modeled as an ARMA process. The latter is due to the impulsive nature of the resulting spectrum from the poles in the Fourier transform of the resulting model. AR methods are based on linear prediction, which is what's used to extrapolate the signal outside of it's known values. As a result, they do not suffer from sidelobes and require no windowing.

Subspace methods decompose the signal into a signal subspace and noise subspace. Exploiting orthogonality between the two subspaces allows a pseudospectrum to be formed where large peaks at narrowband components can appear. These methods work very well in low SNR environments, but are computationally very expensive. They can be grouped into two categories: noise subspace methods and signal subspace methods.

Both categories can be utilized in one of two ways: eigenvalue decomposition of the autocorrelation matrix or singular value decomposition of the data matrix.

Noise subspace methods attempt to solve for 1 or more of the noise subspace eigenvectors. Then, the orthogonality between the noise subspace and the signal subspace produces zeros in the denominator of the resulting spectrum estimates, resulting in large values or spikes at true signal components. The number of discrete sinusoids, or the rank of the signal subspace, must be determined/estimated, or known ahead of time.

Signal subspace methods attempt to discard the noise subspace prior to spectral estimation, improving the SNR. A reduced rank autocorrelation matrix is formed with only the eigenvectors determined to belong to the signal subspace (again, a model order problem), and the reduced rank matrix is used in any one of the other methods.

Now, I'll try to quickly cover your list:

PSD using Burg method: The Burg method leverages the Levinson recursion slightly differently than the Yule-Walker method, in that it estimates the reflection coefficients by minimizing the average of the forward and backward linear prediction error. This results in a harmonic mean of the partial correlation coefficients of the forward and backward linear prediction error. It produces very high resolution estimates, like all autoregressive methods, because it uses linear prediction to extrapolate the signal outside of its known data record. This effectively removes all sidelobe phenomena. It is superior to the YW method for short data records, and also removes the tradeoff between utilizing the biased and unbiased autocorrelation estimates, as the weighting factors divide out. One disadvantage is that it can exhibit spectral line splitting. In addition, it suffers from the same problems all AR methods have. That is, low to moderate SNR severely degrades the performance, as it is no longer properly modeled by an AR process, but rather an ARMA process. ARMA methods are rarely used as they generally result in a nonlinear set of equations with respect to the moving average parameters.
PSD using covariance method: The covariance method is a special case of the least-squares method, whereby the windowed portion of the linear prediction errors is discarded. This has superior performance to the Burg method, but unlike the YW method, the matrix inverse to be solved for is not Hermitian Toeplitz in general, but rather the product of two Toeplitz matrices. Therefore, the Levinson recursion cannot be used to solve for the coefficients. In addition, the filter generated by this method is not guaranteed to be stable. However, for spectral estimation this is a good thing, resulting in very large peaks for sinusoidal content.
PSD using periodogram: This is one of the worst estimators, and is a special case of Welch's method with a single segment, rectangular or triangular windowing (depending on which autocorrelation estimate is used, biased or unbiased), and no overlap. However, it's one of the "cheapest" computationally speaking. The resulting variance can be quite high.
PSD using modified covariance method: This improves on both the covariance method and the Burg method. It can be compared to the Burg method, whereby the Burg method only minimizes the average forward/backward linear prediction error with respect to the reflection coefficient, the MC method minimizes it with respect to ALL of the AR coefficients. In addition, it does not suffer from spectral line splitting, and provides much less distortion than the previously listed methods. In addition, while it does not guarantee a stable IIR filter, it's lattice filter realization is stable. It is more computationally demanding than the other two methods as well.
PSD using Welch's method: Welch's method improves upon the periodogram by addressing the lack of the ensemble averaging which is present in the true PSD formula. It generalizes Barlett's method by using overlap and windowing to provide more PSD "samples" for the pseudo-ensemble average. It can be a cheap, effective method depending on the application. However, if you have a situation with closely spaced sinusoids, AR methods may be better suited. However, it does not require estimating the model order like AR methods, so if little is known about your spectrum a priori, it can be an excellent starting point.
PSD using Yule-Walker AR method: This is a special case of the least squares method where the complete error residuals are utilized. This results in diminished performance compared to the covariance methods, but may be efficiently solved using the Levinson recursion. It's also known as the autocorrelation method.
Spectrogram using short-time Fourier transform: Now you're crossing into a different domain. This is used for time-varying spectra. That is, one whose spectrum changes with time. This opens up a whole other can of worms, and there are just as many methods as you have listed for time-frequency analysis. This is certainly the cheapest, which is why its so frequently used.
Spectral estimation: This is not a method, but a blanket term for the rest of your post. Sometimes the Periodogram is referred to as the "sample spectrum" or the "Schuster Periodogram", the former of which may be what you're referring to.

If your interested, you may also look into subspace methods such as MUSIC and Pisarenko Harmonic Decomposition. These decompose the signal into signal and noise subspace, and exploits the orthogonality between the noise subspace and the signal subspace eigenvectors to produce a pseudospectrum. Much like the AR methods, you may not get a "true" PSD estimate, in that power most likely is not conserved, and the amplitudes between spectral components is relative. However, it all depends on your application.

Cheers


Does a signal and its autocorrelation ever have the same Fourier spectra? My understanding is that Fourier power spectra give the frequencies and corresponding intensities/amplitudes contained within a signal. In certain techniques, such as molecular dynamics, it is common to take the autocorrelation of a time-series before calculating the Fourier spectrum.

I have read that in molecular dynamics this practice can be due to the appearance of autocorrelations in quantum mechanics, and I have also read that taking the autocorrelation removes the dependence on the initial conditions.

a) Does the autocorrelation contain the same frequency components as the original signal?

b) Are there any conditions under which the Fourier power spectrum of the signal and it's autocorrelation become mathematically identical? eg. infinite signal length, perhaps.

c) What are the factors that cause these two spectra to differ?

d) Are there any other reasons to favor using the autocorrelation in place of the original signal when computing the power spectrum?

Thanks.

(a). Yes, a signal and its autocorrelation have the same frequencies. If $X(f)$ is the Fourier transform of the signal and $A(f)$ the Fourier transform of its autocorrelation, then $A(f) = |X(f)|^2$. So, in general, a signal and its autocorrelation don’t have the same Fourier transform, but there is at least one exceptional case:  $x(t) = W \operatorname{sinc}(Wt)$ whose Fourier transform is $X(f) = \operatorname{rect}\left(\frac fW\right)$ and so $A(f) = |X(f)|^2 = X(f)$. Those who wish to generalize this notion to signals whose Fourier spectra are the sum of non-overlapping $\operatorname{rect}$ functions should write out the details for themselves.

All of the above is for finite-energy signals.

If one is thinking of periodic finite-power signals with Fourier series $x(t) = \sum_{n=-\infty}^\infty c_n \exp(j2\pi nf_0t)$, then the periodic autocorrelation function has Fourier series $\sum_{n=-\infty}^\infty |c_n|^2 \exp(j2\pi nf_0t)$. Fourier series are represented by impulse trains in the frequency domain:
\begin{align}
X(f) &=  \sum_{n=-\infty}^\infty c_n \delta (f-nf_0),\\
A(f) &= \sum_{n=-\infty}^\infty |c_n|^2 \delta (f-nf_0)
\end{align}
and if we wish to continue to use the formula $A(f) = |X(f)|^2$ with the spectra of these periodic signals, then we must assume that
$$c_n \delta(f- nf_0)c_m^* \delta^*(f- mf_0) = \begin{cases}
|c_n|^2 \delta(f- nf_0), & \text{if}~m=n,\\
0, &\text{if}~m\neq n.
\end{cases}$$
Note that we are assuming that $|\delta(f-nf_0)|^2 = \delta(f-nf_0)$ and $\delta(f- nf_0) \delta^*(f- mf_0)=0$ for $m\neq n$, but if we are willing to swallow this bald and unconvincing narrative, then everything is hunky-dory.  Once again, we see that the signal and its autocorrelation have the same frequencies but that their Fourier transforms are generally not equal unless it so happens that $|c_n|^2 = c_n$ for all $n$. Note that this means that for each $n$, $c_n$ must equal $1$ or $0$. Thus, the only real-valued periodic signals whose Fourier transforms equal the Fourier transform of their autocorrelation function are of the form
$a_0 + \sum_{n=1}^\infty a_n 2\cos(2\pi f_0 t),~~ a_i \in \{0,1\}, i = 0, 1, 2, \ldots$


Why bit rate is same in a lower frequency spectrum and higher spectrum? If there is a spectrum between 0 Mhz - 1 Mhz, bandwidth would be 1 Mhz here and for example with baseband we can send 2000 bps. What I can't understand is if spectrum is 5-6 Mhz instead of 0 - 1 Mhz, although the spectrum is the same 1 Mhz, shouldn't we be able to send more data because we're sending it with much higher frequencies? Why is bit rate fixated on spectrum instead frequencies itself? What am I missing here?

If you are a bass and sing A 220 Hz and B 247 Hz at 120 beats per minute, you can sing data at a certain rate.

If you are a soprano and sing two quarter notes A 880 Hz and B 988 Hz at the same BPM, your data rate isn't any higher.

If you sing more notes (higher bandwidth) you could communicate a more complex score (e.g. a higher data rate).  So bandwidth (range of notes) matters, not being a tenor or alto-soprano (spectral frequency).

Same at RF frequencies.


Can you average out focal length from different camera calibration sessions? If I perform a number of camera calibration routines at a fixed resolution, with various sizes of calibration patterns (like checkerboards), at various distances and get various values for focal lengths,  could I just average them out? Would this just give a more "robust" focal length? Is it a good idea?

No, it is not a good idea.  Typically, you want to limit the volume of space relative to the camera in which you want to do your measurements, and you want to use a single calibration target appropriate for that volume.  In theory, you can calibrate with any number of different patterns.  In practice, however, most implementations assume that the pattern is the same in all calibration images.

In any case, if you insist on using multiple patterns, you still want to jointly optimize the parameters, not just to average them after the fact.  And you would most likely have to implement that yourself.


Discrete time Signal Transformation I have a discrete-time system defined by  $y[ n]=x[5n+2]$. And an  input $x_1[n]=\delta[n-5]$.

Is the output  $y_1[n]=\delta[25n-5+2]$?

The input to the system $x_1[n]=\delta[n-5]$ is zero everywhere and is only one when the argument of delta is zero (at $n=5$). The output $$y_1[n]=x_1[5n+2]=\delta[(5n+2)-5]$$ would be zero everywhere too, except when the argument of $\delta$ is zero: $$5n+2-5=0$$ which requires $n=\frac{3}{5}$. But since $n$ only has integer values, it can never be achieved. Hence, the output of the system is zero always.


Need of $\tt abs()$ method when plotting a power spectral density for a given dataset I am a newbie in signal processing and would like to know the significance of using the abs() function and squaring the values received as an output of fft() function of Python's scipy.fftpack library, used when trying to plot a power spectral density for a dataset. I have found that many of code examples to plot a power spectral density do use an abs() and then square the values obtained thereafter.

Can anyone please provide me a reason for doing so?
Can't we just directly plot the values obtained from fft() function in Python's scipy.fftpack library?

Here is the code I have written till now to plot a power spectral density by referring some of the code examples,

import scipy.io as sio
import numpy as np
Import matplotlib.pyplot as plt
import pandas as pd

df = pd.read_csv("denoised.csv")
data = df.values

from scipy.fftpack import fft,fftfreq
dft= fft(data)
PSD = np.abs(dft) ** 2

Complex sines, or cisoids, $e^{iw}$ are fundamental functions for the study of linear  systems, even if the latter are real. And the FFT is a fast algorithm for a discretized version of the continuous Fourier transform.

So when you apply an FFT, it provides you with a bunch of coefficients $c_k$, which can be positive, negative, or complex. We are often interested in the energy carried by these coefficients, defined as their squared modulus, which can be computed  "complex number times conjugate equals square of modulus", or $c_k\times \overline{c_k} =|c_k|^2 $. This works with real numbers too. So instead of squaring up, you could as well multiply the FFT by its complex conjugate, pointwise.

Calling it "absolute value" is a bit of an abuse. This term is generally understood for real numbers only, and "modulus" would be a better term. However, abs is a very common term in this situation.

Note that sometimes, people are interested in the argument of the complex coefficients, and instead of 2D plots use the third dimension to "unwrap" the phase part, sometimes called phasor/polar representations:


Is my concept of receive diversity right? EGC,MRC and SC Is my  concept of receive diversity right?There are three combining technology belong to receive diversity.

1.Equal gain combining(EGC)

$\sum e^{-j\theta_k}(\sqrt{P}|h_k| e^{j\theta_k}x[n]+w[n])$

so the EGC coefficient is $e^{-j\theta_k}$ just because we assume the channel phase is $e^{j\theta_k}$,so it seems that this method is hope the signal can have  same phase direction and sum it.

2.Selection combining(SC)

we have a receiver with multiple antenna,so theoretically we have multiple SNR in the receiver,now we select the highest SNR as real SNR. It seems that this method like Ostrich Syndrome,that is ,we don't do any other process to the signal,we just choose the highest SNR and ignore the other weaker SNR.

3.Maximal ration combining(MRC)

we use the Cauchy–Schwarz inequality to realize the MRC coefficient is the conjugate of the channel,that is

$\sum |h_k|e^{-j\theta_k}(\sqrt{P}|h_k| e^{j\theta_k}x[n]+w[n])$

So theoretically,the SNR performance order is

MRC>EGC>SC,because MRC not only do some process in the phase,but also amplitude,however,EGC just do some process in the phase.we all know that the signal power has a high relation with amplitude,so the performance of MRC will be the highest in these three method.

Is my concept and thinking right?

Sounds pretty solid to me.

It's not necessarily MRC>EGC>SC, you'd need to write $\geq$ technically, since for the special case where all $|h_k|$ are equal you have MRC=EGC and for the special case of a single branch (antenna), they are all equal.


regarding plots in visual dsp ++ 5.0 For my source code witten in c ,in visual dsp++ 5.0 for a forward error correcting code on adsp 21062 sharc family platform. i want to plot a graph between different data sets. when i tried to use the plot window, it is showing the data sets which are in assembly.so how to plot a graph between the input and output in the tool?please help me out.

phani tej

export the data from memory and process it in MATLAB, Octace of whatever your favorite analysis tool is. VDSP has only very limited data analysis capabilities


Noisy Frequency Estimation using Hilbert Transform I'm looking for a way to estimate instantenous frequency with low time-variation. The following MATLAB function (save as test_fest.m) does this using the Hilbert transform on a single 10 Hz sine tone. Sample rate f_s is 1000 Hz.

function test_fest()
%TEST_FEST()  Test Frequency Estimation using Various Techniques.

fS = 1000;
f = [10]';                                  % Sine Wave Frequency
nFFT = fS;
x = sin(linspace(zeros(length(f),1), f*2*pi, fS));
X = hilbert(x,nFFT);

fE = diff(unwrap(angle(X))) * fS / (2*pi); % Frequency Estimate
eE = mean((fE/f-1).^2); % MSQ Error

plot(fE,'.');
title(sprintf('Hilbert Transform Instantaneous Frequency Estimation of Sine Wave f=%.2f Hz, f_S=%g Hz, Mean-Square Error:%.2g', f, fS, eE))

Output is shown here:

I, however, get a greater error, 2.3e-5, than I expected. Is there something wrong with my approach?

Digital Hilbert transform is tricky to do since an ideal Hilbert transformer requires an FIR filter with an impulse response from -infinity to +infinity. Obviously that can't be done, so numerical problems are unavoidable. They tend to be worst at low frequencies so my guess is you'd see a smaller error at, say, 100 Hz. Andrew Duncan's paper "The Analytic Impulse" is thorough analysis of what the digital Hilbert transform can and can't do. See http://www.andrewduncan.ws/air/index.html

Some alternative methods of measuring frequency have been discussed here: Is there an algorithm for finding a frequency without DFT or FFT?


Biquad cascading frequency response I am implementing EQ lines in my filter. I need to do biquad with slope 24 db/oct, 36/oct and bigger so I would like to cascade several biquads. I am stuck with calculating frequency response of two biquads - now I am able to calculate magnitude and phase and draw one curve (peak filter).
Should I multiply this response by itself? Thank you in advance for the answers.

The application of a filter to a signal in the time domain is defined as convolution.  Convolution in the time domain is mathematically equivalent to multiplication in the frequency domain. So if we filter a signal, and then filter the output of the first filter, this is equivalent to multiplying the signal by both filters in the frequency domain. Therefor, this is equivalent to filtering with a third single filter whose frequency response is equal to the product of the first two filters’ frequency responses. Because frequency responses are complex numbers, the resulting magnitude is equal to the product of the two input magnitudes and the resulting phase is the sum of the input phases.  This can be easily extrapolated for more than two filters as well.


Hypothesis testing with binary pattern Let's say you have a sequence $\{x_i\}_{i = 1}^N$ of ones and zeros. You know that $P(x_i = 1) = \frac{1}{3}$.
You want to test two hypotheses
$H_0$: The sequence is $iid$ with $P(x_i = 1) \ \forall \ i$
$H_1$: The sequence is full of $001$ sub-sequences, but has noise. For example, the sequence might start off as $001001$, then have a period of random noise (where $P(x_i = 1) = \frac{1}{3}$).
After the period of noise, the sequence matches the $001001$ pattern again. Importantly, the $1$ values might land on different values than when it started.
That is, we are not testing whether $mod(i, 3) = 0 \implies x_i = 1$. We are merely testing whether the "normal" state of this sequence is batches of $001001001$ but those batches could start at any value.
Another way of putting this is that we allow for "extra" bits to be inserted.
Ideally, I would like to find a test statistic that enables me to reject the null in this scenario. One possible solution, suggested in the comments below, is to count the number of $001$ sequences, and compare to the number of $001$ sequences you would find under the null.
This is an interesting solution. I will have to think more about how one might construct a confidence interval or assymptotics for this statistic.
I appreciate your comments.
Update: After thinking about the scenario that we are working with more, I think I can safely assume that only $0$ values are inserted. This is a much simpler problem.
$H_0$: The sequence is $iid$ with $P(x_i = 1) = \sum_{i = 1}^N x_i$. We guess it's Bernoulli distributed according to the number of $1$s we observe in the data.
$H_1$: The sequence has lots of $001001001$ segments, but extra $0$ values get inserted at random into the sequence (which creates offsets as above).
I don't think this materially changes the approach too much, but means we no longer need to worry about the distribution of the "Noise" sequences.

Given what you explained, you can count the number of times the sequence is
001, 100 or 010. Then, $\hat{p}$ is defined as that count divided by the number of sequences, $N$.
If we define the null as the hypothesis that the proportion of good sequences  is 1.0,  then the test is
$H_0: p = 1.0$ versus $H_1: p \neq 1.0$
As far as testing the null, ( I'm assuming N is large here so that the normal
approximation to the binomial holds )  you can calculate the test statistic
as $Z =  \frac{ \hat{p} - p_{0}}{\sqrt{\frac{\bar{p}(1-\bar{p})}{N}}}$
and use $p_{0} = 1 $ and $\bar{p} = 0.5$. Then, reject the null if $|Z| > Z_{\frac{\alpha}{2}}$ where $\alpha$ is the critical value of the test.
The CI  would be $\hat{p} ~ \pm Z_{\frac{\alpha}{2}} \times
{\sqrt{\frac{\bar{p}(1-\bar{p})}{N}}}$.
If a digit in one of the three digit sequences is is missing or there's an extra, and you know this, then you can just count that sequence as a bad sequence.


Why only odd harmonics after non-linear amplification? Why do I get only odd harmonics appearing? This not only happens at saturation of the transfer curve, but also in the linear portion of the curve.

I make a input signal that has too high power that it will saturate the transfer curve. The transfer curve below is for Amplitude modulation, I also have a transfer curve of the output phase change for each input signal level.

This is basically how I do it

My input signal, x(t), is the sum of two sine waves

For each point on the transfer curve above I apply Amplitude and Phase modulation by using this array

modulating_phasor[] = (power_out / power_in) * exp(i*phase_deviation)

Then I look at my input signal and lookup where its value matches the input power of the transfer curve and use that location in the modulating_phasor array

output_signal[] = x[t] * modulating_phasor[location]

This makes the output signal complex array so I take the absolute value and for parts that should be negative, I multiply that location on by -1 by finding the correct location on input signal.

I cant upload my code... but you get the idea of what I am doing right?

Now this is my output and no matter what my input level on my signal is i.e. whether I am in the non-linear part or linear part, I always get odd harmonics. I've also tried a different transfer curve. Only get odd again.

An ideal amplifier would have a transfer characteristic of $f(x)=Ax$: the input signal comes out amplified and otherwise undistorted. A real amplifier will deviate from this and go into saturation. We could model it by a polynomial $f(x) = \sum_{n} a_n x^n$. Now, what we would still expect is that the amplifier treats positive and negative values the same way, i.e., it is symmetric in the sense that $f(-x) = f(x)$. This leads to the requirement $a_n = 0$ for all even $n$.

Now, one can show that raising a sinusoidal with frequency $f_0$ to the $n$-th power generates harmonics up to order $n$, however, preserving their oddity:

$n = 2$ generates frequencies 0 and $2f_0$
$n = 3$ generates frequencies $f_0$ and $3f_0$
$n = 4$ generates frequencies 0, $2f_0$ and $4f_0$
$n = 5$ generates frequencies $f_0$, $3f_0$ and $5f_0$

and so on. It's a consequence of the binomial theorem and quite easy to show (see, e.g., here).

In short: A transfer function that is odd symmetric will generate only odd harmonics. This is why you are seeing only those.


QPSK demodulator on GNURadio continuously loses phase lock I am designing a QPSK demodulator to be used for receiving Meteor M2 (LEO satellite) signals (QPSK at 72Kbaud). However, I am finding a lot of problems in my way.
First of all, this is the flowchart I came with:

First, I filtered the signal, applied it a $\frac{3}{16}$ fractional downsampling and fed it to a Polyphase Clock Sync block (8sps input, 2sps output and 2pi/100 bandwidth, 0.5 TED Gain and 0.7 Damping factor and 0.15 Deviation). Then, since I thought that for lower passes multipath would be a thing, I added a CMA equalizer after the symbol sync block. Lastly, I placed a Costas Loop just before the soft decoder.
When the elevation of the satellite is above ~50°, I obtain a pretty clear constellation with no decoding errors. Nevertheless, the interesting challenge of designing such demodulator is for it to work at lower elevations. It is very important that it does not lose lock at any moment of the pass, since even the shortest lock loss will produce a black line on the final image.
I am having a lot of trouble with the signal lock at low elevations. I assume that this is caused by a "dirty" signal at the input of the Polyphase Clock Sync block.
• Does my flowchart look good? This is the first digital receiver I design in GNURadio and have very little experience in this field.
• Is there any way (other than trial and error) to optimise Symbol Sync block parameters?
• What technique (filtering, equalization, etc.) may I use in order to solve the locking problems?
Thank you in advance!
Edit: the input signal to the low pass filter has already no doppler shift.
Edit 2: can a PLL Carrier Tracking block help in this situation? The signal seems to be frequency stable, but I do not know whether such block can help avoiding phase locks at the symbol sync stage. Can placing the Costas Loop before symbol sync help?
Edit 3: after some refinement, this is the look of the receiver

Thank you very much @gotchi85 por the help

From what I understand, the modulation is a QPSK at 72k symbol.
So the main lobe is 144 kHz wide.
First point: Shaping filter clarification
Do you know if there is a shaping filter used by the transmitter?
Usually it is an RRC filter with a roll off of a given value.
If there is one in the transmiter, your RRC filter in the receiver should match roll off value
Second point: The low pass filter
The main lobe is 144kHz wide (because of QPSK-72ksps). The low pass filter you use is 90 kHz (2 x 45 kHz, for negative + postive spectrum side) wide.
With a 0.5 alpha roll off used at the transmitter, its bandwith would be 72 x (1+0.5) = 108 kHz <=> 54 x 2. With a 45 kHz, you cut before a possible RRC inside the chain (45 vs 54). Expand the low pass filter bandwith so that it does not interfere with the RRC filter in the receiver of Symbol block.
Try 60 instead of 45 Khz
Third point: AGC

AGC placement:
Is the AGC placement intended?
What happens if the signal is too low at the input of the first filter?
Do you have an external AGC before entering GNU radio?

AGC rate:
I have a small experience with GNU radio, so don't know what is behind the rate or other AGC blocks...
The important thing is that the AGC does not track on the varation of the symbols but on the average signal level.
Try maybe to decrease that rate so the reponse is slower.

There are probably other points to tune, these are main I see. You could earn several fractions of dB in those first places.
If there is an improve, you should see it on the constellation directly (if not too noisy)
The best approach is to first make the design and verify it performances in lab conditions and not directly on the satellite.


Appropriate channel bandwidth for my BPSK modulated Psuedo-Random code on direct conversion transceiver My hardware is an analog device SDR+Xilinx FPGA (sdr = fmcomms4-ebz). The direct conversion transceiver (ad9364) has a tunable channel bandwidth of 200Khz - 56Mhz. I am employing BPSK modulation on a PRN code of 100K bits in length (repeats itself).My Tx Sampling rate is 6 Msps, thus the highest frequency present in my transmitted signal (barring harmonics) is 6Mhz.

My Questions:

1) According to Nyquist criterion the minimum bandwidth should be (Fs/2) where Fs is the sampling frequency. However, the main lobe of the BPSK spectrum is 2*Fs wide. Is my minimum Tx bandwidth 3Mhz or 6Mhz ?

2)On the Receiver side if I am oversampling by a factor of 2 (ie Fs_Rx = 2*Fs_Tx) should I adjust my Receiver's channel bandwidth to the sampling rate of the receiver or the transmitted signal? ie should my Rx channel bandwidth be equal to my transmitters channel bandwidth or double ?

3) Beyond the appearance of higher order harmonics and inefficient use of bandwidth, are their negative effects to using an arbitrarily large bandwidth for both my Tx and Rx nodes? (effects like ISI,RF interference)

EDIT 1:

Pulse shape currently using root raised cosine filter, however I am also testing unfiltered data i.e. rectangular pulses. How does this effect the bandwidth (i assume the root raised cosine would allow for lower bandwidth)

Inquiring about both (sorry I should have specified). My current understanding is that the required bandwidth at baseband is half that of passband. Is that correct?

Thanks !

It seems from your question that you are still lacking a few concepts and techniques that you'll have to master before succeding. I'll just give a few pointers:

Get a good book on SDR. I recommend "Software Receiver Design" by Johnson et al, and "Software Radio: A Modern Approach to Radio Engineering" by J. Reed.
The passband bandwidth is double the baseband bandwidth. However, your radio works with quadrature signals and complex sampling, which muddles that concept a little, since the baseband quadrature signal's spectrum is no longer symmetrical.
Do not use rectangular pulses in an RF application.
The actual bandwidth depends on the pulse shape (and in the case of root raised cosine, its rolloff factor). If your sampling rate is $f_s$, then your baseband bandwidth goes from $-f_s/2$ to $f_s/2$.
The receiver's sampling rate is independent of the transmitter's, because the receiver sees an analog signal. The only thing you need to worry about is meeting Nyquist and having enough samples for your synchronizers to work.
You may need to worry about radiating large amounts of energy on a wide band. You should restrict your transmissions to one of the ISM bands and in general comply with all the relevant regulations in your country.
Regarding transmit bandwidth: Let's say you want to transmit at baud (pulse) rate $R_p$, using (root) raised cosine pulses with rolloff factor $\beta$. Then, the baseband signal bandwidth is $$B = \frac{R_p(\beta+1)}{2}.$$ The transmitter sampling rate needs to be larger than $2B$, but for practical reasons it is usually a good idea to oversample at least by a factor of 2. Also, you'll need to select a frequency that is compatible with your hardware. As an example, if I wanted to transmit at $R_p = 100,000$ pulses per second with $\beta=0.75$, I'd calculate $B=87.5$ kHz, and then oversampling by 4 I'd end up with $f_s = 350$ kHz. I'd configure the hardware for the smallest supported sampling frequency that is larger than 350,000.


Software-based sound card OFDM I'm trying to implement an sound card-based OFDM modem. A paper says:

  ...A serial to parallel converter is applied and the IFFT operation is
  performed on the parallel complex data.

Say I have a stream $S$ of complex numbers, then after serial-to-parallel (S/P) conversion, it becomes $S_1$ and $S_2$. But as far as I'm concerned, OFDM requires the sub-carriers to be orthogonal to each other. So can I just perform IFFT twice directly on the streams without defining the frequencies of sub-carriers? I'm quite puzzled. Any help on this will be greatly appreciated. Thanks.

The serial/parallel converter shown in most OFDM transmitter block diagrams is required in a hardware implementation but not very meaningful for software implementation. Because in software there is no IFFT block that actually works with parallel inputs, it's all serial processing in reality.

Let $S(k)$ be your input stream of complex numbers and $N$ be the number of subcarriers, i.e. the IFFT size. Then you have to calculate the IFFT of every $N$ samples in $S(k)$, i.e. of the vectors $[S(iN),\ldots,S((i+1)N - 1)]$, with symbol index $i=0,1,2,\ldots$. The result is the OFDM modulated time domain signal. If I understand your example correctly you have to apply the IFFT to both $S_1$ and $S_2$, yes.

The output of the IFFT operation is also complex, in general. Unless you use two audio channels followed by an Inphase/quadrature mixer that shifts the signal to a carrier frequency the vectors defined above have to exhibit a complex conjugate symmetry in order to obtain a real-valued time domain signal. (So-called discrete multitone (DMT))

The allocation of data to subcarriers happens through the IFFT algorithm you're using. The FFTW library, for example, allocates $S(0)$ to frequency zero and $S(N-1)$ to "frequency" $-N/2$. What physical frequency this corresponds to in turn depends on the sampling frequency of the digital-to-analog converter you're using. Orthogonality is inherent to the (I)FFT algorithm, you don't have to apply any further processing to obtain it.


Help with Hankel transform filters I found a filter for Hankel Transform . As far as i understand this is effectively used as an approximation for $J_0(x)$ and $J_1(x)$.

But how will the change in input affect these filter terms (e.g. $J_0(r) \rightarrow J_0(2r)$)

References : http://onlinelibrary.wiley.com/doi/10.1046/j.1365-2478.1997.500292.x/full

I think you misunderstand what this filter is doing. It does not compute approximations of the Bessel functions $J_0(x)$ and $J_1(x)$, but it computes an approximation of the Hankel transform of some input function. The function is expressed as a weighted sum of Bessel functions, and the result of the transform are the weights. The digital filtering approach is an efficient alternative to solving the transform integral by standard methods of numerical integration.


How does the Meyr-Oerder Timing Algorithm work? I`ve came across a Paper for this Algorihm for estimating the Symbol-Timing from Meyr and Oerder, which I want to understand. The Algorithm uses blocks from the samples and calculates the Timing Estimation for the difference between the Transmitter-Clock and the Receivers clock. By squaring the blocks we are getting a peak in the spectrum which is our symbol timing.

My questions are:

Why are we getting the symbol timing by squaring the signal block?
What does the Fourier-transformation do to the squared signal-block?

For PAM signals, the squaring creates a spectral component at $\dfrac{1}{T}$.  Consider the following trigonometric identity to understand why this is the case:

$$\cos^2\left(2\pi \dfrac{1}{2T}\right) = \dfrac{1}{2}+\dfrac{1}{2}\cos\left(2\pi \dfrac{1}{T}\right)$$

The phase of that spectral component at frequency $\dfrac{1}{T}$ is a measure of the evenness (or oddness) of the symbol relative to the reference time $t=0$ of the fourier transformed squared samples.  A phase of zero corresponds to the symbol being centered at that reference time.

The Fourier transform is used to find the phase of the $\dfrac{1}{T}$ spectral component, which is an unfiltered, scaled estimate of the symbol timing offset.


Relation between FIR filter's transition width and phase delay How does the transition width of a FIR filter relate to phase delay at the output of the filter?

I have been trying to find information on the subject for days.

In this answer I gave two well-known heuristic design formulas relating the filter order of (linear-phase) FIR low-pass filters to the transition bandwidth. According to both formulas, the transition bandwidth is inversely proportional to the filter order. Since the phase delay of a linear-phase FIR filter equals $M/2$, where $M$ is the filter order, the phase delay $\tau_p$ and the transition bandwidth $\Delta\omega$ are also inversely proportional to each other:

$$\Delta\omega\propto \frac{1}{\tau_p}$$


Removing the 'red' artifact from my JPEG I have a JPEG that originally came from a IP Camera. I have converted the JPEG to a bytes and have uploaded it to my server. Upon receiving it on my server I save it to disk.

When I view the image it looks 'grainy' and 'reddish'.

When I perform a contrast stretch I get this:

When I perform a Hist-Equalization I get this:

N.B.
I have also tried a filter just on the Red Channel:

Image<Bgr, byte> CurrentFrame = new Image<Bgr, byte>((Bitmap)picOriginal.Image.Clone());
CvInvoke.cvSmooth(CurrentFrame[2], CurrentFrame[2], SMOOTH_TYPE.CV_GAUSSIAN, 3, 3, 1.5, 0);

Is there another approach so I could make the image look better, 'balanced', if that makes sense?

Did you try to change the white balance of your image ? Apparently there's a function (balanceWhite) in OpenCV that tries to correct the white balance, you might want to give it a go.


Proving orthogonality How to prove a cosine signal with irrational period and another cosine signal with a rational period are orthogonal?

For example, $\cos(4t)$ and $\cos(4\pi t)$ sum is aperiodic. Hence it is said that these two signals are orthogonal over infinite limits? But how to prove this? Normal integration of product of the two signals cannot be used because we get $\cos(\infty)$ which is unknown.

Sinusoids are power signals, and in the case of power signals, cross-correlation between two signals $x(t)$ and $y(t)$ is defined by

$$R_{xy}(t)=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}x(\tau)y(t+\tau)d\tau\tag{1}$$

If $x(t)=\cos(\omega_1 t)$ and $y(t)=\cos(\omega_2 t)$ then for $\omega_1\neq\omega_2$ the cross-correlation at $t=0$ is

$$\begin{align}R_{xy}(0)&=\lim_{T\rightarrow\infty}\frac{1}{2T}\int_{-T}^{T}\cos(\omega_1 t)\cos(\omega_2 t)dt\\
&=\lim_{T\rightarrow\infty}\frac{1}{2T}\left[\frac{\sin((\omega_1-\omega_2)T)}{\omega_1-\omega_2}+\frac{\sin((\omega_1+\omega_2)T)}{\omega_1+\omega_2}\right]\\&=0
\end{align}$$

Note that the definition (1) avoids all convergence problems that you would encounter otherwise.


Intuition behind commutativity of convolution in LTI systems Why is convolution commutative, as it seems to treat two signals in a different way in an LTI system?

If you imagine $y[n] = x[n] \star h[n]$ with $x[n]$ being an input signal and $h[n]$ being the impulse response of an LTI system A, how does it make sense that LTI system B with input $h[n]$ and impulse response $x[n]$ generates the exact same output $y[n]$ ?

In a discrete-time system such as the one that you have, the number $y[n_0]$ (here $n_0$ is a fixed integer) is a sum of the form $$\sum_{k=-\infty}^\infty h[k]x[n_0-k]$$ which can be re-arranged via a change of variables (replace $k$ by $n_0-\ell$) to $$\sum_{\ell=-\infty}^\infty h[n_0-\ell]x[\ell].$$ So, the commutativity of the convolution is trivial.  The issue is the interpretation that you put on it.  As Laurent Duval/s answer points out, the systems A and B are not equivalent in any sense of the term.  If the signal $x$ were replaced by a different signal $\hat{x}$, then system A would have output $\hat{y} = h \star \hat{x}$, but you wouldn't get the same output $\hat{y}$ if system B were excited by $h$; the impulse response of system B continues to be $x$, and system B thus has output $x \star \hat{x} =  \hat{x}\star x \neq h \star \hat{x}$.


Plotting frequency spectrum for a signal with 20MHz bandwidth and sampling frequency of 12.5GHz I am new to Matlab and I am trying to plot the frequency spectrum of an OFDM signal with bandwidth of $20 \text{ MHz}$ and DAC sampling rate of $12.5 \text{ GHz}$. The number of samples is $400,000$. Plotting the signal I am only able to see the Nyquist bandwidth but not the $20 \text{ MHz}$.

I would appreciate any help. Link to the data file.

Here is the frequency spectrum plotted via Matlab :

Below are my line of codes:

 BW = 20MHz;
 f_DAC = 12.5GHz;
 Nfft = 32;
 Sig_fft = abs(fftshift(fft(Data_signal)/length(Data_signal))).^2;
 Freq = ((0: length(Sig_fft)-1)-length(Sig_fft)/2)*(f_DAC/length(Sig_fft));
 figure();
 plot(Freq, 10*log10(Sig_fft/1e-4));
 title('Signal Spectrum');

Your fundamental limitation is the (comparatively) short window length of $M=400,000$ samples, at such a prohibitively high sampling frequency of $12.5$ GHz. This will roughly yield a few hundred kHz of effective FFT resolution per bin of the FFT spectral samples.

Using the formula for the spectral resolution of various windows we can deduce the following resolution, for example, for a standard Hamming window as

$$\Delta w = 8\pi / M$$ where $M$ is the window's length, that is $M=400,000$, which yields $$\Delta w = 8\pi / 400,000 = 6.2832 \times 10^{-5}$$ radians.

This resolution is determined by the window length, and its frequency equivalent under the sampling rate of $f_s = 12.5 \times 10^9$ Hz is

$$ \Delta f = \frac{\Delta w}{2\pi} f_s = \frac{6.2832 \times 10^{-5}}{2\pi} 12.5 \times 10^9 = 125,000 Hz$$ i.e., the effective frequency resolution of your window is $\Delta f = 125$ kHz. Now given your bandwdith of $B=20$ MHz, this equates to a number of $20 \times 10^6 / 125\times 10^3 = 160$ samples.

On the other hand the apparent count of the number of FFT samples (apparent resolution) inside the $20$ MHz region yields something about $k = 4\times 10^5 \times 12.5\times 10^9 / 20\times 10^6 = 640$ samples, which however does not carry more information than those inferred from $160$ subsamples.

Hence a bandwidth of about $20$ MHz would yield about less than 200 true-samples of a spectral region inside $400,000$ many apparent-FFT samples. This is quite narrow. You would either use a much longer observation interval or use a much reduced sampling frequency as @MBaz commented.


Basic FIR Filtering with audio buffer in C++ I am willing to implement a basic FIR filter using C++ and JUCE.

I wrote this simple algorithm using the information I found in the DSPGuide, but I only seem to get a variation in volume for all the frequencies, instead of a frequency filtering.

How can I correct it to effectively apply frequency filtering ?

if (channel == 0) {

    float result = 0; // initialisation to 0 of the result

    for (int sample = 0; sample < buffer.getNumSamples(); sample++) { // for each sample

        // moving the samples in the delay line
        for (int i = 0; i < 8; i++) { // delayLine being 9 values long
            delayLine[i + 1] = delayLine[i];
        }
        delayLine[0] = reader[sample];

        result = 0;
        for (int i = 0; i < 9; i++) { // for each tap
            result = result + delayLine[i] * filterTaps[i]; // multiply
        }
        writeback[sample] = result; // output
    }
}

Here are the coefficients used, which are supposed to produce a Low-Pass Filter in the audible range :

filterTaps[0] = 0.002385;
    filterTaps[1] = 0.011910;
    filterTaps[2] = 0.026352;
    filterTaps[3] = 0.039825;
    filterTaps[4] = 0.045351;
    filterTaps[5] = 0.039825;
    filterTaps[6] = 0.026352;
    filterTaps[7] = 0.011910;
    filterTaps[8] = 0.002385;

I have generated them using this website http://www.arc.id.au/FilterDesign.html
Can't remember the cut-off I chose, but it was between 1000 and 5000 kHz with a 44100 Hz sampling rate. These coefficients are just an example, as my initial goal was to make an anti-aliasing filter. Since the filter was not working properly, probably because of the algorithm, I chose another set of coefficients in the audible range to ease the debug process).

Thank you.

Your implementation of the delay line is flawed.  It's just copying the second last sample over each element.

Try this instead:

    for (int i = 7; i >=0; i--) { // delayLine being 9 values long
        delayLine[i + 1] = delayLine[i];
    }

Previous answer below

Your filter seems to be a lossy low pass filter --- even the passband has -15dB attenuation.

Try using these coefficients instead:

0.01156111 0.05773286 0.12773940 0.19304879 0.21983567 0.19304879 0.12773940 0.05773286 0.01156111

They are obtained by normalizing the sum of the coefficients you provided.

Your filter has this response:

The normalized coefficients give this response:

R Code Below

 #30506
filterTaps <- c(0.002385,
0.011910,
0.026352,
0.039825,
0.045351,
0.039825,
0.026352,
0.011910,
0.002385)

filterTaps2 <- filterTaps / sum(filterTaps)

freqz(filterTaps)

freqz(filterTaps2)

getNumSamples <- 1000

reader <- rnorm(getNumSamples,0,1)
writeback <- 0*reader

delayLine <- 0*filterTaps
result <- 0; # initialisation to 0 of the result
for (sample in seq(1,getNumSamples))
{ # for each sample

    # moving the samples in the delay line
    for (i in seq(8,1,-1))
    { # delayLine being 9 values long
        delayLine[i + 1] = delayLine[i];
    }
    delayLine[1] = reader[sample];

    result = 0;
    for (i in seq(1,9))
    { # for each tap
          result = result + delayLine[i] * filterTaps[i]; # multiply
    }
    writeback[sample] <- result; # output
}


Can a wavelet transform give time dependent phase of sinusoids in signal I have a signal which contains sinusoidal components that oscillate at different frequencies. I think the phase of the sinusoids is changing with time. I could do a Fourier transform on small chunks of the signal and the phase of the sinusoids from that, but is there a better way to do this with a wavelet transform? I'm new to the concept of wavelets. From what I understand, it's good at doing what an FFT does but with time resolution. Could you get the phase of a sinusoid as a function of time from a wavelet transform ?

It is indeed possible (up to a period), with complex "continuous" wavelet transforms, provide your mother wavelet ooscillate sufficiently. You can look at the scalogram modulus and its phase at the appropriate scales, as illustrated from Matlab: Wavelet coherence:

Alternatively, you can convert your signal into an analytic one (using the Hilbert transform) and perform a real wavelet analysis.

Finally, you might dig into other linear (short-term Fourier transform) or bilinear time-frequency (Wgner-Ville, Choi-Williams) or time-scale tools.


What do I do with the result of an FFT? So I am working on a project using digital signal processing, which I am not familiar with at all, and a lot of this is all new to me. I am using SciPy\NumPy to play around with code to understand the Fourier Transform.  Basically, I have been able to read in a wave file and perform an FFT on the data using the following code:

samplerate, data = wavfile.read(sys.argv[1])
fftres = fft(data)

for i in range(1, len(fftres)):
    print i, ",", fftres[i], "\n"

My results look something like this:

(51.0478396084-421.235956694j)
(0.0324635664874+372.476395255j)
(-198.202964067+219.743659247j)

I'm not sure what I'm looking at to be honest.  I need to find the frequencies in my file so I can apply a bandpass filter on the audio signal.  From my output, how do I actually get the frequencies?

I would recommend doing some proper reading of books/tutorials etc on the Fourier transform and the Discrete Fourier Transform (DFT).  You will also want to look at filters and probably convolution for the bandpass filter.  But regardless, here's some info that should help with your current state:

Use matplotlib to plot the FT that you've calculated:

fs = 44100 # Or whatever the actual sample rate is (Hz)
f = linspace (fs,len(data), endpoint=False)
mpl.plot (f, abs (fftres))
title ('Magnitude spectrum of the signal')
xlabel ('Frequency (Hz)')

Another function you will want to look at is fftshift, which makes it simple to 'roll' your spectrum around so that DC is in the centre.

The frequencies that are sampled by the DFT are between $-f_s/2 < f < f_s/2$ for sampling frequency $f_s$.  And the FFT produces the same number of samples as the number of data points provided to it.  The FFT function returns the spectrum so that the DC (constant/average) value is first in the array.  Use the function fftshift to adjust the array if you want it ordered  $-f_s/2 < f < f_s/2$.

I will now repeat: Fourier transforms and filtering are a huge field, and while it will definitely be interesting and useful to play with the functions for a bit like you are doing, you really need to read up on what is going on if you want to actually do anything useful.

Aside: why are you starting at index 1 in your loop, rather than 0?


How to get the volume level from PCM audio data? I have a decoded PCM audio data obtained from the Web Audio API.

I'd like to generate a chart that represents the volume level (perceived by humans) over time.

volume|
 level|    ******
      |   *      *                           **
      |  *        *                         *  **
      |**          *      ***              *
      |             ** * *   *            *
      +---------------*-*-----************------+--- time
      0                                        30s
          audio is             and quiet
          loud here            here

Here is my current approach, given I have an audio composed of 6 channels.

To calculate the volume level at time T:

For every channel, calculate the RMS value of 200ms window (44100 * 0.2 = 8820 samples, in my case). Time T will be in the middle of the window (ignore the edges).
Calculate the average of the 6 values obtained in step #1.

I get something like this:

Is this a reasonable approach to get the volume level?

Would translating RMS to dB (step #3: Calculate 20*Log_10(step #2)) produce a more accurate chart?

What would be some appropriate methods for smoothing / noise removal?

There are two separate issues behind your question: the duration of the measurement and the algorithm used to do the actual metering. In the audio jargon, the term used for the perceived level measurement is loudness metering.

Time duration

The first issue you need to tackle is the length of the window used to do the metering. Usually, in audio they can be divided in three classes:

Momentary, which gives a loudness level that varies quickly over time (usually a few hundreds of milliseconds)
Short term, which winows of a few seconds
Integrated, which gives a loudness level on the whole audio file from beginning to end

The actual values for the lengths of the windows varies from the specification, but this EBU tech note specifies a duration of 400 ms for momentary, 3 s for short-term and the whole duration of the file for integrated, of course. Moreover, if you do live tracking of the metering, there is the issue of the refresh rate of the measurement as well, but you may want to set this aside in your case.
Finally, it may be useful to note that a rectangular window is used, contrarily to other cases in audio (short-term fourier transform for instance).

Loudness metering

To accurately represent the perceived level, you will need something more complex than a simple RMS measurement. Various methods have been developed over time, from the simple A-weighting mentioned in the comments (which is just filtering before computing the RMS value), to the LUFS (Loudness Units relative to Full Scale) measurement specified in the ITU-R BS.1770-4 recommendation. The latter is a complex algorithm, but it provides methods for measuring perceived loudness in a multi-channel audio file, which is your case.

Regarding your last questions:

Translating the measurement to dB is a good idea to produce a more understandable chart, and most audio products that implement loudness metering give the value in dB.
The windowing method already has inherent smoothing, because giving one value for a window of, say, 400ms time is already a kind of low-pass filtering. Regarding noise removal, most modern algorithms implement gating, i.e. an algorithm that removes all noise below a given level, which means that long periods of silence in the audio file will not lower the value measured when doing integrated measurements. It is then possible to compare the loudness from different sources (between a 2-hours movie and a 3-minutes rock song for instance).

Finally, please note (as it has been highlighted in the comments) that the value will be given in dB relative to full-scale, so it will not reflect the level of the actual audio source, except if you know the calibration of the digital-to-analog converter used for the initial recording.

There are a lot of audio plugins that implement metering, and one of the most recent ones is the Fabfilter Pro-L. You may not need such complexity for your measurement, but it can bee a good idea to take a look at what is done in the industry, the way the value is given and what a good loudness curve looks like.

You can find more comprehensive and extensive sources on loudness metering on the web, and especially on web sites of pro audio manufacturers, for instance TC Electronic or Waves


Relation between order and stability in IIR filter What is the effect of increasing order of IIR filter? does it also effects stability? If we have an IIR filter of order 6, now we change its order to 7, will there be any change in its stability?

Stability depends on the pole locations, so one cannot generally say that increasing the filter order does or doesn't affect stability. It depends on how you increase the filter order. What are the additional coefficients of the denominator? (Because the denominator polynomial of the transfer function determines the pole locations).

Another question are the finite word length effects. If you quantize the filter coefficients, e.g., by using fixed point arithmetic (but of course also by using floating point), then higher filter orders tend to be more problematic than lower filter orders, because the chance that poles get very close to (or even outside) the unit circle through quantization becomes larger. But again, you cannot say anything in general unless you specify the filter's transfer function, the filter structure, and the type of quantization.


Removing white pixels from a b/w image in MATLAB? Consider this MATLAB code:

mg=imread(........bmp');
level=graythresh(mg);
x=im2bw(mg,level);
figure,imshow(x);
BW2 = bwperim(~x);
figure,imshow(BW2);

operating on this image (NOTE has white margin, see here):

The result was (original link here):

How can I remove the white pixels from the second image and then obtain the first image with everything but the white pixels removed from the second image.

In other words, how can I remove the boundaries from the second image and then get the original image without these boundaries?

So...you want to remove the fat white borders from around the black region??

You could try search for the first and last row,column that are 0 in the second image, using (e.g.) find:

% assuming BW2 is 2D, one channel.
% the reason the first has a BW2' and the second is BW
% is because of Matlab being column-major when it does `find`:
% have to make sure we actually get the first/last indices.
[col1,row1] = find(BW2'==0,1,'first'); % top-left corner
[row2,col2] = find(BW2==0,1,'last'); % bottom-right corner
BW3 = BW(row1:row2, col1:col2);
imshow(BW3)

This does assume that the top-left and bottom-right corner of the image proper won't be part of the outline - this is a very low chance.

If you want to be sure however, you could instead look for the first row that is not all 255 and first col that is not all 255, and similarly for the last row/col (ie locate margin by looking for entire rows/cols being 255) and crop like that:

cols = sum(abs(255-BW2),1); % work out difference from 255
rows = sum(abs(255-BW2),2);
col1 = find(cols,1,'first');
col2 = find(cols,1,'last');
row1 = find(rows,1,'first');
row2 = find(rows,1,'last');
BW3 = BW(row1:row2,col1:col2);

(Note this will fail on your jpeg because of the compression it uses -- the border is not uniquely 255, there are a few 251's in there. But it will work on your code if you don't save BW out as a jpeg first and then read it back in. If you're going to save images, use a lossless format like TIF or PNG).


How to estimate ALL time-delays across MULTIPLE signals? Imagine there are 100 recordings of the same signal or pattern, but those recordings are not properly aligned with each other in time. In other words, each sample has some unknown time delay in relation to each of its sister samples. The data might look something like this:

The question is: What's the best method for computing the relative time-delays of those signals so that I might align them?

If there were only two signals, computing time-delay would be simple. The best method I'm aware of is to maximize the cross-correlation between the two samples across different delay values.

But what about computing a large set of delays? I'm sure there has to be a robust method for this general task... I just don't know what it's called. Thanks!

There are papers describing a method called multichannel cross-correlation, that recursively estimates time delay of each signal.

However, If I have not missunderstood you problem you can also select one of your signal to be your reference signal and then repeatedly calculate cross-correlation against that compared to the other signals. that way you can align every signal to your reference point.

EDIT: based on your latest comments:

  So far I've only had a chance to look at the original paper linked (microphones in a sound room). It's very interesting, however it seems to require that you know something about the geometry of the microphones. I had been hoping to find a generalized, system-agnostic method to do multi-channel cross-correlation. Something with which you can just input a set of arbitrary signals or patterns and it outputs the "best fit" delay vector.

What it appears to me that you need then is  blind system identification methods which I can link you some examples. Common methods for this task utilizes cross correlation error minimization but there are some limitations based on your model selection.

Under-Modelled Blind system identification - cross correlation based

Blind system Identification for acoustic source localization - Eigenvalue decomposition method

Blind system identification using sparse leaning method - LMS adaptive algorithm based on cross correlation formulation


3D image sample frequency I am currently frequency filtering a 3D image in matlab, with x,y,z dimensions of different lengths. The matlab filter takes normalised frequency as an argument. What does this mean for an image and how do I calculate it?

Note : the distance between each pixel is 2mm

2 mm period between samples is 500 samples per meter.  That's your sampling frequency.

Matlab's digital filter design functions (is that what you're using?) take a normalized frequency from 0 to 1, where 1 is the Nyquist frequency (1/2 the sampling frequency; the highest spatial frequency you can measure at that sample rate).

So if you wanted to filter at a spatial frequency of 100 cycles per meter (features that cycle every 1 cm), for instance, you would note that the Nyquist frequency is 250 cycles per meter, and enter a value of 100/250 = 0.4 for the normalized frequency.


ARMA model vs. Linear Filter ARMA Model

$$y_t = \sum_{k=1}^p a_{k}y_{t-k} + \sum_{k=0}^q b_{k}\epsilon_{t-k}$$

where $b_0$ is usually 1 and $y_t$ is observable.

General Linear Filter

$$x_t = \sum_{k=1}^p a_{k}x_{t-k} + \sum_{k=0}^q b_{k}y_{t-k}$$

$y_t$ is observable (input), $x_t$ is filter output, and there is a feedback loop via $x_{t-k}$.

ARMA is a model, but from filtering perspective it looks like a filter for $\epsilon_t$, by simply moving $y_t$ to the right side of the equation and $\epsilon_t$ to the left side of the equation, then normalizing by $b_0$ (usually 1 so no effect).

Is a model always a filter at the same time?
Given a linear filter can we directly infer the underlying process model by simply rearranging the terms?

For example given an FIR filter to recover $x_t$ from $y_t = x_t + \epsilon_t$

$$\hat{x_t} = \sum_{k=0}^3 b_{k}y_{t-k}$$

Can we say that the whole process is $AR(3)$ with Gaussian errors given the filter is optimal in MSE sense?

$$y_t = \frac{b_{1}}{-b_0}y_{t-1}+\frac{b_2}{-b_0}y_{t-2}+\frac{b_3}{-b_0}y_{t-3}+\epsilon_t$$

This doesn't look correct.

I think the following information is relevant, but didn't really clarify it for me.

https://www.mathworks.com/help/signal/examples/linear-prediction-and-autoregressive-modeling.html

  Linear prediction and autoregressive modeling are two different problems that can yield the same numerical results. In both cases, the ultimate goal is to determine the parameters of a linear filter. However, the filter used in each problem is different...In the case of linear prediction, the intention is to determine an FIR filter that can optimally predict future samples of an autoregressive process based on a linear combination of past samples. The difference between the actual autoregressive signal and the predicted signal is called the prediction error. Ideally, this error is white noise. For the case of autoregressive modeling, the intention is to determine an all-pole IIR filter, that when excited with white noise produces a signal with the same statistics as the autoregressive process that we are trying to model.

So according to this text, the corresponding model should actually be $AR(4)$

$$y_{t+1} = \sum_{k=0}^3 b_{k}y_{t-k} + \epsilon_t$$

  two different problems that can yield the same numerical results

Are there situations where they wouldn't?

Based on your exploration of the topics I can make the following explanation, mostly following the approach in Statistical Digital Signal Processing and Modeling, HAYES.

Linear Time Invariant (LTI) filtering is associated with the following LCCDE (Linear Constant Coefficient Difference Equation) with zero initial conditions and assuming a causal solution; i.e., interested in the solution for $n \geq 0$ in general.

$$ \sum_{k=0}^{N} a_k y[n-k] = \sum_{k=0}^{M} b_k x[n-k]  \tag{1} $$

where $a_k$ and $b_k$ are the coefficients of the LCCDE (mostly $a_0 = 1$), and $x[n]$ is the input, while $y[n]$ is the output of the filter (aka solution of LCCDE).

The filter is associated with an impulse response $h[n]$ (coefficients of the LTI filter) which is the solution of Eq(1) when $x[n] = \delta[n]$; $x[0] = 1$ and $x[n] = 0$ for all $n> 0$.

$$ x[n] \longrightarrow \boxed{ h[n] } \longrightarrow y[n] $$

Impulse response can be used to compute the solution $y[n]$ of the LCCDE in Eq(1) for a given input $x[n]$, through what's called a convolution sum :

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

Z-transform of the impulse response is the following and is useful as an algebraic transform representation of the LTI filter system.

$$ H(z) = \frac{ \sum_{k=0}^{M} b_k z^{-k} }{\sum_{k=0}^{M} a_k z^{-k}} \tag{3}$$

From which the output $y[n]$ can also be expressed as
$$ Y(z) = H(z) X(z) \tag{4} $$

which is nothing but Z-transform applied to Eq(2).

ARMA-(p,q) modeling :

Consider an LTI filter with LCCDE coefficients $a_k$ and $b_k$, as given in Eq(1), then if the input of this filter is a white-noise $v[n]$ with unit variance, then the output of this filter is called (by definition) Auto-Regressive Moving Average process with orders (p,q), abb as ARMA-(p,q), the associated LCCDE is written like

$$ \sum_{k=0}^{p} a_k x[n-k] = \sum_{k=0}^{q} b_k v[n-k]  \tag{5} $$

Note that with the filtering convention in Eq(1), $x[n]$ was the input and $y[n]$ was the output, but withing the ARMA-(p,q) convention in Eq(5), $v[n]$ is the input and $x[n]$ (the ARMA process) is the output. The AR order is $p$ (associated with the lefthand side, recursive part, of the LCCDE) and the MA order is $q$ (associated with the right hand side of the LCCDE).

Now given some WSS RP $x[n]$, then you can model this process as an ARMA-(p,q) by finding the set of coefficents $a_k$ and $b_k$ according to Eq(5), which would minimize the mean square error between the modeling ARMA-(p,q) filter output $\hat{x}[n]$ and the given process $x[n]$ to be modeled. The solution for $a_k$ are generated by Modified Yule-Walker Equations and $b_k$ are found in a more complex manner following MA modeling techniques etc.

Linear Prediction :

Given a WSS random process $x[n]$, we want to estimate its current value $x[n]$ by  using its $p$ past values  $x[n-1],x[n-2],...,x[n-p]$ via a linear combination :

$$ \hat{x}[n] = \sum_{k=1}^{p} c_k x[n-k] $$

I have used these $c_k$ as the linear prediction coefficients; it can be shown that these coefficients are related to all-pole signal model Prony coefficients, or an ARMA-(p,0) model coefficients as $c_k = -a_k$.

Rewriting Eq(5) with $q=0$,and $a_0 = 1$, as

$$ x[n] = - \sum_{k=1}^p a_k x[n-k] + b_0 v[n] $$

with a suitable $b_0$, reveals that, if the predicted process $x[n]$ was generated via an ARMA-(p,0) model using Eq(5), then the linear prediction coefficients $c_k$ which  minimize the prediction error will be $c_k= -a_k$ of Eq(5).


Can realizations of a filtered Gaussian white noise process be represented as a Fourier transform? Suppose we have a noise process $V(t)$ which is the result of passing Gaussian white noise through a filter with frequency response function $H(\omega)$.
Can we represent realizations of this process as a Fourier transform
$$V(t) = \int \frac{d\omega}{2\pi} M(\omega) e^{i \phi(\omega)}$$
where $M(\omega)$ and $\phi(\omega)$ are random variables?
If so, what are the statistics of $M(\omega)$ and $\phi(\omega)$?

I don't think it makes sense to talk about the Fourier transform of a (realization of a) random process, because it generally doesn't exist if no further restrictions are imposed. This problem is also discussed here.

What can be done is analyze the statistics of the discrete Fourier transform (DFT) and the discrete-time Fourier transform (DTFT) of a windowed version of a discrete-time stationary random process. This has been done in this paper.


The Wiener Filter for the Time Derivative of the Signal Assume the observation of the original signal $s(t)$ is $x(t) = s(t) +n(t)$, where the signal and noise are independent.
Then we need to design the Wiener filter $g(t)$ to estimate $\frac{d s(t)}{dt}$ which is the derivative of the original signal $s(t)$.
The design of the Wiener filter is from the MMSE: $E[e]$. But how to formulate the error function here? $e(t) = \frac{d s(t)}{dt} - x(t)$ or $e(t) = \frac{d s(t)}{dt} - \frac{d x(t)}{dt}$?
I think this can be generalized to any linear operators like estimating $y = Hs$ (deconv, deblur etc.). But I have no idea the begining step. Could anyone give me some hints? Thanks!

The key point here is the derivation / integration is a Linear Operator.
Since the Wiener Filter is also a linear operator it makes things easier.
Write the problem using the Wiener Filter Model:

Neglecting the integrator we know the Wiener Filter for $ \hat{H} \left( w \right) $ (See Wikipedia - Wiener Deconvolution):
$$ G(f) = \frac{1}{\hat{H}(f)} \left[ \frac{ 1 }{ 1 + 1/(|\hat{H}(f)|^2 \mathrm{SNR}(f))} \right] $$
By defining $ H \left( w \right) = \frac{1}{j w} \hat{H} \left( w \right) $ as the fourier operation of integration (This is not exactly accurate, see The Integration Property of the Fourier Transform)
$$ G(f) = \frac{1}{{H}(f)} \left[ \frac{ 1 }{ 1 + 1/(|H(f)|^2 \mathrm{SNR}(f))} \right] $$
Basically, we defined the integration as a filter and used the known building blocks.
The interpretation is, with respect to the non integral case, is apply the derivative operator on the estimator and have some regularization on the SNR estimation as we estimate a derivative.
Resources

How Is the Formula for the Wiener Deconvolution Derived?
Question on Wiener Filtering.
1D Deconvolution with Gaussian Kernel (MATLAB).


How do I calculate the autocorrelation of a train pulse having only the pdf? I have a pulse train of a fixed period, but a random amplitude following a Gaussian distribution $(\mu=0,\sigma=1)$. I know how to work with a random period or delay, but how do I proceed to calculate the autocorrelation with a random amplitude?

Suppose that your pulse train is denoted $X(t) = Y\cdot g(t)$ where
where $g(t)$ is the basic pulse train and $Y$ is the random
amplitude of the pulse train.  Then, $X(t)$ is not a
wide-sense-stationary process and its autocorrelation function is
$$R_X(t_1,t_2) = E[X(t_1)X(t_2)] = E[Y^2]\cdot g(t_1)g(t_2).$$
If $Y$ has mean $\mu$ and variance $\sigma_2$, then
$E[Y^2] = \mu^2+\sigma^2$.  If $g(t)$ is a function of period $T$,
that is, $g(t) = g(t+T)$ for all $t$, then $R_X(t_1,t_2)$ is
periodic in both arguments:
\begin{align}
R_X(t_1+T,t_2) &= (\mu^2+\sigma^2)g(t_1+T)g(t_2)
= (\mu^2+\sigma^2)g(t_1)g(t_2) = R_X(t_1,t_2)\\
R_X(t_1,t_2+T) &= (\mu^2+\sigma^2)g(t_1)g(t_2+T)
= (\mu^2+\sigma^2)g(t_1)g(t_2) = R_X(t_1,t_2)
\end{align}
I am quite sure that this is not what you hoped to hear and so you
must revise your question appropriately, giving more details and
describing more clearly what it is you want to find out.


Relationship between sensors and materials in terms of electromagnetic spectrum The following text is from Digital Image Processing(Gonzalez and Woods)[3ed, Ch 2.2,p46]

  If a sensor can detect energy radiated by a band of electromagnetic spectrum, we can image events in that band. But the wavelength of electromagnetic wave required to “see” an object must be same size or smaller than the object.For example, a water molecule has diameter $10^{-10}$ m. To study water molecules, a source capable of emitting in far ultraviolet or soft X-ray region is required.

For the water molecule example, why do the authors refer to a source emitting ? As an analogy, for X-rays, the bones absorb waves in X-ray region and this renders them opaque on the film. Therefore, to image water molecules, a source which emits far ultraviolet/soft X-rays is needed. Passing this through water would somehow let us image water molecules. Is that what they mean ?

But the wavelength of electromagnetic wave required to “see” an object must be same size or smaller than the object.

Essentially, the waves are unaffected (at least not enough to be detected) by objects that are smaller than the wavelength. Since we look at the distortion in the waveform to image objects, no distortion would mean an object is undetectable.

Imagine you are walking down a rocky path. Small rocks and pebbles would not change your gait at all, so they would be undetectable to you. Large stones, however, would force you to change your step, making them detectable. If you were much smaller and took shorter steps (or a shorter wavelength), then these small objects would become detectable to you.

So, to answer your questions:

  For the water molecule example, why do the authors refer to a source
  emitting ?

Because it is the wavelength of the radiation source which matters for imaging.

  Therefore, to image water molecules, a source which emits far
  ultraviolet/soft X-rays is needed. Passing this through water would
  somehow let us image water molecules. Is that what they mean ?

At that wavelength or smaller, yes. Operating at that wavelength, with the proper detection equipment, should allow you to image water molecules and larger.


How to sketch the following discrete time signal $\sum_{k=0}^{n} \delta[n-k]$? I tried to solve it like in the picture but It feels so wrong though.

Your signal is:
$$
x[n]  = \sum_{k=0}^n \delta[n - k]
$$
and you're saying you want to sketch it.

Based on your comment, I'm going to assume that
$$x[-1] = \left . x[n] \right  |_{n=-1} =  \left . (\delta[n] + \delta[n+1]) \right  |_{n=-1} = 0 + 1 = 1$$
and that
$$x[+1] = \left . x[n] \right  |_{n=+1} =  \left . \left(\delta[n] + \delta[n-1]\right) \right  |_{n=+1} = 0 + 1 = 1$$

I believe that means
$$
x[n] = 1, ~~ \forall n \in \mathbb{Z}
$$
because
$$
x[n]   = \cdots + \delta[n - (n-2)] + \delta[n - (n-1)] + \delta[n - n] + \cdots
$$
if summing up or
$$
x[n]   = \cdots + \delta[n - (n+2)] + \delta[n - (n+1)] + \delta[n - n] + \cdots
$$
if summing down.

It's a little odd having the upper limit of the summation depend on the current time and, as I said in the comment, having the summation go down to the upper limit if the upper limit is below the lower one.


Common methods of transmitting digital data through an extremely noisy audio signal? I am developing software for transmitting data from a computer to a mobile device using audio. Specifically, from the computer's speaker to the device's microphone. (I am developing both the sender and the receiver software.) This is one channel, and the sample rate is a standard 44.1Khz. It needs to be VERY resilient to noise - it will be used in environments with street noise, people talking, background music, etc., which may be 10dB louder than the signal itself.

Knowing almost nothing about DSP, thus far I have apparently implemented frequency-division multiplexed on-off keying, without realizing it until now :) Basically each of 96 frequencies corresponds to a bit, and on the receiving end, if the frequency is present, the bit is on. I then use some custom error correction code to whittle those 96 bits down to 8 very reliable bits.

My system works very well and is very robust, but unfortunately has an EXTREMELY low data rate (~80 bits/sec). I have some ideas about how to improve bandwidth, like using more than 96 frequencies, and using the phase component to encode additional data. Also it appears I should look into frequency-shift keying and phase-shift keying.

So - my question is simply this:

What are the options for encoding data in a noisy broadcast audio signal, and what kind of bandwidth and error rates can be expected with each technique?

The scheme you are using is called On/Off Keying.  It is not terribly efficient, but it is simple and gets the job done.

When you say that the signal is 10 dB below the noise floor I suspect what you mean is that if you add up all the signal energy and all of the noise from 0.3 - 14 kHz the signal is 10 dB weaker, but that the signal uses a much narrower frequency range, making the signal stronger at very specific frequencies.  If it wasn't, you would not be getting reliable data transmission without doing long synchronous averaging, and it sounds like you aren't.

There are a number of different ways you could do this.  The state of the art would be to use an OFDM signal in the 17-21 kHz band.  A much simpler but not as efficient way to go (but still much more efficient than your OOK) would be to use a QPSK signal in the 17-21 kHz band.  If you have enough transmission power you could modify it to 8-PSK or 16-QAM for higher throughput.  Since it sounds like the transmission power is more of a limitation than the bandwidth (you can get a lot more data through 4 kHz of bandwidth than 80 bits/s), I would add a forward error correction (FEC) code to the data stream.

With QPSK you can get (ignoring the potential reduction in throughput from the error correction code) 1 bit/s/Hz, so you would have an upper limit of 4 kbits/s.  With 8-PSK it would be 6 kbits/s, and with 16-QAM it is 8 kbits/s.  The price of going with the higher data rates is that you are more likely to have errors.

The data rate for OFDM is more complicated.  Given your level of DSP knowledge I would go with QPSK, and then modify from there if needed.  It sounds like you would be thrilled with anything even close to 4 kbits/s, so that should be plenty good enough.


why derivative and logarithm operators amplify the noise effects? It is said, do not use derivative operator in a the control system. because it will amplify the noise effect. Can someone explain it and give some mathematical reference about it please?
Also I heard, we shouldn't use logarithm operator too, for the same reason. How come?
Thanks in advance.

The derivative operator is linear and time-invariant and can be described by its frequency response

$$X(j\omega)=j\omega\tag{1}$$

From (1) it is obvious that its gain increases linearly with frequency, and high frequencies (where the SNR is usually bad) are greatly amplified (as already pointed out in endolith's comment).

The logarithm is a memoryless nonlinearity and it compresses the input signal. This operation is used in many signal processing algorithms, such as feature extraction for speech recognition. In these algorithms this operation is essential to the performance and it definitely does not deteriorate the signal. What is important, however, is to deal with the singularity of the logarithm as its argument approaches zero because: $$\lim_{x\rightarrow 0+}\log(x)=-\infty$$
If nothing is done about it, the logarithm will produce very high negative peaks during periods where the signal is almost zero (i.e. when there's only noise). This problem is easily handled by defining a threshold, either on the value of the logarithm or on the input signal. Note that the logarithm can only be applied to non-negative signals (usually magnitudes or squared magnitudes of frequency domain data).


Help Understanding Time of Flight Phase Sanitization Algorithm I'm implementing the technique found in this paper: http://front.sjtu.edu.cn/~jinyh/gradCourse/Paper/Wireless/SpotFi%20Decimeter%20Level%20Localization%20Using%20WiFi.pdf

to perform localization.

I've implemented 95% of the technique, but this one short algorithm is giving me problems.

The explanation in pseudocode is very simple but the results I'm obtaining are wildly different than reported in the paper, and I don't think this is just a case of misreporting results.

Notice that to minimize the function in step 1, rho will always go to 0! Which is clearly not right.

An example of this algorithm being applied in the paper is given here, with the first figure being the unwrapped phase prior to running Algorithm 1 and the second one being the resulting phase.

I can share my code so far for Algorithm 1 but I really think that I'm just missing something pretty essential in the description....

I think it's just a simple line fitting, assume that phase = at + b
then the phase without offset is phase_without_offset = phase - at. I used ployfit function in matlab to implement the algorithm.The following picture is my result.

I guess you are a student from SJTU, the same school with me, i'm also implementing the spotfi, can you tell me your email or qq, i want to talk with you! my email is mrshao at sjtu.edu.cn(please replace at with @).


Is this input function time varying? I'd like to understand when an input signal is time varying or time invariant.
Specifically, I have this function:

$y(n) = x(n) + 0.1 + (0.2)^n + x(n-3)$

I thought that this is time varying because of the $n$ in the exponent. Would anyone mind to clear this up for me? Is this time invariant or time varying?

Signals are (almost) always time variant. Otherwise they would be very boring :-). I believe your question is about whether "systems" are time invariant or not. A system basically is a "prescription" of how to create an output signal from an input signal. The equation you show is a system equation of how to create output y[n] from input x[n].

The idea of time invariance is the following: The behavior of the system doesn't depend on absolute time, i.e. the relationship between input and output is the same, regardless of whether it's Tuesday or Friday. You can find a more mathematical description here http://en.wikipedia.org/wiki/Time-invariant_system

The easiest invariance test is the following:

pick any input signal, typically $\delta [n]$ will do fine.
Calculate the output $y{_{0}}[n]$ for that intput
now delay the input signal by some amount.
Calculate the output $y{_{1}}[n]$ for the delayed input
If $y{_{1}}[n]$ is NOT simply a delayed version of  $y{_{0}}[n]$ then the system time variant

Another (equivalent) test is the following

For time invariant system, the system equations only dependency on time must be as linear index into the input signal

Either one will work.


Find cut-off frequency of a filter given its transfer function I have a high pass filter given by the time-domain relationship between input and output:

$x'(t) = x(t) - \alpha \times x(t-1)$.

I would like to get its cut-off frequency. How can I proceed?

First you discretize you equation (note: I assumed that your $x(t-1)$ means a delay of one samples and not 1 second, but you can adapt easily if that's not the case) :

$$ x'[n] = x[n] - \alpha \times x[n-1] $$

Where
$$ x[n] = x(n\times T_s)$$

Then you use the $\mathcal{Z}$ transform :
$$X'[z] = X[z] - \alpha \times X[z]\times z^{-1}$$

And then you get the transfer function in the Z domain :
$$ \frac{X'[z]}{X[z]} = 1 - \alpha \times z^{-1}$$

This is a causal FIR filter, called a 1-zero filter with its zero at $\alpha$. Your cutoff frequency will depend on your sample frequency.

Illustration :

Code for the illustration :

figure,
hold on
for alpha=0.1:0.1:1
    h = freqz([1 0], [1 alpha]);
    f = 0:1/length(h):1-1/length(h);
    semilogy(f, 20*log(abs(h)));
    hold on
    grid on
end
legend('0.1', '0.2', '0.3', '0.4', '0.5', '0.6', '0.7', '0.8', '0.9', '1');
title('normalized frequency response')


Can two different impulse responses give the same frequency response? I have two different impulse response each with different length, is it possible that they have the same frequency response?

No. The impulse response and frequency response of an LTI system are related by the Fourier transform, which is one-to-one.


Transferring data over sound with QAM / QPSK Been working on a FSK Modem in the last few weeks and have gotten it to work with 400bps, very reliable, atleast over a virtual audio cable, but it also works when playing a sound directly into a microphone.
Now I've shifted my focus to other Modulation techniques that could be used, such as M-QAM (16, 64, 256) and was wondering if anyone has any experience in that direction.
Would QAM be viable for the transfer of data over audio? What concerns me is the complexity of the signal detection when it comes to QAM, or even QPSK. I am trying to do this, same as the FSK, in purely Python, with .wav files as the intermediate medium of storing and playing the signals.
Would appreciate any kind of information or tips.

The problem with transmitting data over an acoustic link is the variability of the channel. Things you need to take into account are

Bandwidth, quality & directionality of the microphone
Bandwidth, quality & directionality of the loudspeaker
Specifically the distortion properties of the loudspeaker
Relative placement and distance of speaker and microphone
Acoustic environment: reverb time as a function of frequency, distinct reflections, absorption properties, sympathetic buzzing and rattling, etc.
Background & competing noises: number of sources, spatial distribution, temporal & spectral properties, energy,
All of the above can be stationary or time variant.

As a result the amount of data you can reliably transmit varies by several orders of magnitude with the channel properties alone. In a clean & stable environment with good speaker close to a microphone QAM is definitely possible. In a busy restaurant with a laptop speaker: not so much.


DSP Topic for Newcomers Does any one of you know what topic to teach people related to signal processing.
The audience is mainly for a group of software programmer working on a telecommunication product who have limited knowledge of signal processing.
Thanks!

I think the  most fundamental things that anybody trying to understand DSP should know are:

Nyqvist Sampling Theorem and Effects of Aliasing
Need for DFT and the use of FFT to achieve the same.
Basics of Digital Filters.

You have mentioned the audience as telecommunication engineers. So you could also add a lecture on the  need for modulation,some basic modulation and demodulation techniques and need for interpolation and decimation.


Are IIR filters (and specifically Butterworth filter) causal? I am new to signal processing. I am learning now about filtering and I am trying to implement high pass filter (HPF),and later band-pass filter (BPF). I implement the filters in Python 3.8 with NumPy and SciPy.
For the offline case (filtering a pre-recorded time) I designed a Butterworth high-pass filter, and applied on the signal at the time domain using signal.filtfilt function. However, as the documentation shows, this approach is good only for offline cases and not for real-time, since the filter scans the signal forward and then backwards. I can recover the frequency response from the filter using signal.freqz function.
Now I want that this filter will work also in real-time (so it needs to be causal, and I ask if this possible for this type of filters. Here are my questions:

If I understand correctly, Butterworth filters are IIR (infinite impulse response) filters. Is it true?
Are IIR filter causal? It may be that all IIR filters are not causal and it may be that some are and some are not. What are the cases in which IIR filters are causal?
Are Butterworth filters causal? If not always, in which case they are causal?
If the answer to 3 is true, how to implement a causal Butterworth filter, both at the time domain and the frequency domain? (I can use SciPy functions such as signal.butter, signal.freqz etc.)

Relevant links:

Help designing Butterworth filter
filtfilt: https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.filtfilt.html

Thank you.

Yes, Butterworth are IIR.  The decay from an impulse technically lasts forever.
Yes, all [implementable] IIR are causal.
Yes, because of #1 and #2.
Don't use signal.filtfilt.  Use signal.lfilter.  filtfilt does the same thing as lfilter, except twice, in opposite directions, which changes a causal filter into a zero-phase filter.

However, as the documentation for both of those functions suggests, for most practical use, you should actually use the SOS variants instead:

The function sosfilt (and filter design using output='sos') should be
preferred over lfilter for most filtering tasks, as second-order sections
have fewer numerical problems.

To summarize the functions:

lfilter: Causal single-stage filtering (low orders only)
filtfilt: Zero-phase single-stage filtering (low orders only)
sosfilt: Causal second-order sections filtering
sosfiltfilt: Zero-phase second-order sections filtering

Really the only reason to use lfilter or filtfilt is if your coefficients are already in b, a format, if you're implementing something out of a textbook, etc.
If you're designing the filter yourself, then just use SOS form, which reduces numerical error.
Note that the filtfilt functions apply the filter twice, so it will have double the order of the original filter.


Designing a Chirp in the Frequency Domain I am trying to create a chirp signal in the frequency domain, by simply setting the magnitudes of each frequency bin to one, and then setting the phase delay function such that the resulting time delay increases for each frequency bin. I am trying to calculate the length of the resulting chirp from the phase delay function.
If the time delay for each frequency, $td(f)$, equals the period multiplied by the number of cycles introduced by the phase delay, and $\phi(f)$ is the phase delay in radians as a function of frequency, then
$$
td(f) = \frac{\phi(f)}{2\pi}\frac{1}{f}
$$
And if we define $\phi(f)$ as,
$$
\phi(f) = 2\pi nf^2
$$
Then it follows that the time delay, $td(f)$, is given by,
$$
td(f) = \frac{2\pi nf^2}{2\pi}\frac{1}{f} = nf
$$
And it seems obvious to me that the maximum time delay would occur when f is at its maximum $\frac{fs}{2}$, where $fs$ is the sampling frequency.
However, when I try and synthesise the signal (in Matlab),
% Create frequency axis (0 to nyquist)
fs = 44100;
fAx = 0:1:(fs/2 - 1);

% Define n so td = nf = 0.25s, and define Phases
n = 0.25/(fs/2);
phi = -2*pi*n*fAx^2;

% Calculate the maximum time delay
td = n*(fs/2);

% Calculate the phase of each frequency up to nyquist
chirp_fft = exp(1i*phases);

% Create Hermitian symmetric signal and ifft
chirp_fft = [chirp_fft(1:end-1), 0, flip(conj(chirp_fft(2:end-1)))];
chirp = ifft(chirp_fft);

% Create time axis and plot chirp
tAx = 0:1/fs:(length(chirp) - 1)/fs;
plot(tAx, chirp)

So by my calculation, the chirp should have a maximum delay of 0.25 seconds, however, when I plot it:
It is 0.5s long. When I try this for other values, the resulting chirp is consistently twice as long as my calculation says it should be. Practically this is not a problem, but I'd really like to understand where I've gone wrong! Where is my missing factor of two?

Instead of phase delay $td(f)$, use group delay $\tau(f)$:
$$\tau(f) = -\frac{d}{d\omega}\phi(f),\tag{1}$$
calculated as the negative of the derivative of the phase $\phi(f)$ with respect to the angular frequency $\omega$, which is defined by:
$$\omega = 2\pi f\quad\Leftrightarrow\quad f = \frac{\omega}{2\pi}.\tag{2}$$
Using your phase shift $\phi(f)$:
$$\phi(f) = 2\pi n f^2 = \frac{n\omega^2}{2\pi}, \tag{3}$$
we get from Eq. 1:
$$\tau(f) = -\frac{d}{d\omega}\frac{n\omega^2}{2\pi} = \frac{2n\omega}{2\pi} = 2nf.\tag{4}$$
This should explain your empirical results, considering that here $\tau(f) = 2\,td(f)$.


Alamouti code after or before IFFT block in OFDM We are trying to implement Alamouti STBC on underlying ofdm signal. We were thinking as to where to write alamouti coder block.

I thought we could write it after the cyclic prefixer block as we get a symbol after the IFFT and cyclic prefixer. we could take two such symbols,send both of them at one time instant and then conjugate them and send them on the next time instant.

But then we saw a few block diagrams where alamouti coder is used  before IFFT block and after ofdm symbol mapper. when I searched why that's the case, there were reasons that said for OFDM to be effective, symbols should be assigned to subcarriers in frequency domain to mitigate multipath effects and not in time domain. The points that I don't understand are:

What does it mean by saying symbols in ofdm  are assigned in frequency domain on individual subcarriers? I can only see a symbol at the output of a ofdm mapper which is a complex number which comes from the constellation diagram depending on the modulation technique used. The output of the IFFT is also a complex number. So  how does an IFFT operation done on complex symbols mean that they are assigned in frequency. we could as well use a FFT instead of a IFFT.
What is the reason for alamouti code to be implemented before IFFT? Can it also be implemented the way I thought after cyclic prefixer block?

The basic concept of OFDM is to divide a high-bitrate datastream into $N$ low-bitrate datastreams and to multiplex these low-bitrate datastreams in frequency. That is, every datastream is assigned to a distinct frequency band (so-called subcarrier) that does not interfere with the other frequency bands. Orthogonality allows the frequency bands to be packed close together which increases the spectral efficiency of OFDM. As it turns out the frequency multiplex can be implemented very efficiently using the IFFT algorithm that transforms a frequency-domain signal into a time-domain signal. Every $N$ complex QAM symbols are one input vector to the IFFT which produces a vector of $N$ time-domain samples that are then transmitted sequentially over the channel. The construction of the transmit signal in frequency domain has several advantages:

Fast fading channels (often the case for multipath mobile channels) might significantly attenuate some frequencies. OFDM enables adaption to such a channel by switching off some subcarriers or decreasing the modulation order on the respective subcarriers.
Usage of cyclic prefix as guard interval makes channel equalization easy to implement.
Dynamic resource allocation is possible (see LTE uplink)

Your suggestion to take $N$ QAM symbols, apply the FFT and then transmit the FFT output corresponds to transmitting the spectral analysis of some time-domain signal one bin per symbol interval. This is totally different from applying the IFFT and I don't see any advantage in doing so.

  What is the reason for alamouti code to be implemented before IFFT? Can it also be implemented the way I thought after cyclic prefixer block?

Let's assume we apply the Alamouti space time block code (STBC) after the cyclic prefix. Let $x_k$ be the discrete-time output of the cylcic prefix block. Then the coded signal vector is given by
$$
\mathbf{x} =
\left(
\begin{matrix}
  x_{k+1} & -x_k^*\\
  x_k     & x_{k+1}^*
\end{matrix}
\right)
$$
If we assume a perfect channel with channel matrix $\mathbf H = (1\quad1)$ and if we ignore noise, the received signal is
$$
\mathbf r = \mathbf H \mathbf x =
\left(
\begin{matrix}
  x_{k+1} + x_k & -x_k^* + x_{k+1}^*
\end{matrix}
\right)
$$
It is clear that we first have to decode $\mathbf r$ before we can demodulate it with the FFT. For doing so, the channel matrix has to be known. In general, the channel matrix is unknown at the receiver and thus has to be estimated. This channel estimation is done in frequency domain in OFDM systems, i. e. after the FFT block. As a consequence the STBC decoder has to be placed after the FFT as can be seen from the following figure:

(Figure is taken from Kim, J et al.: Receiver Designs for Alamouti Coded OFDM
Systems in Fast Fading Channels)

This in turn means that the STBC encoder has to be placed before the IFFT.


Two-sided impulse response and its inverse in real life Dear StackExchange gurus,

I would like to pose this question from a practical standpoint rather than a theoretical one, although perhaps some context in the latter might come useful.

Real-life devices (say Digital-to-Analog converters for example) have a roll off in their frequency response due to analog front-ends. To find said response (and their inverse) one can use spectral methods (T(f) = out(f)/in(f)), or time-domain methods that will extract the impulse response (Wiener-Hopf equations, LMS, etc).

Now, in time-domain methods we can chose our solution to have the dominant coefficient at the center tap of the filter, or we can choose it to be the first tap (i.e. single sided and with zero delay). Both of these give similar results in some occasions, but there are times when the single sided inverse does not invert the system as well as the two-sided does.

Can someone give me some insight on why this happens? Which of the two approaches is better for equalization (regardless of implementation complexity)?

The analog front end will always contain some sort of low pass filter. This can be chosen to be either linear phase or minimum phase. Linear phase preserves the phase response of the original signal but the time domain response is non-causal, i.e. the impulse response is symmetric around the center tap.

Minimum phase has a one sided impulse response, with the strongest tap up front, but it incurs some phase distortion. Ideally we want both, but we have to pick one.

The same consideration implies to the inverse: the inverse of a linear phase filter is linear phase and the inverse of a minimum filter is minimum phase. However there is an added complexity: not all filters are invertible: If the magnitude of the original transfer function is very small or zero at some frequencies, the inverse is not well defined. You get "divide by zero" or "divide by noise" problems.