FFT command line application? Are there any command line programs for Windows, preferably free and stand-alone, which can report the peak/strongest frequency within a given range of frequencies?
I need something like this to automate finding the frequency of a calibration signal which slowly drifts.

Estimating Multiple Frequencies

The Goertzel algorithm is more efficient than the DFT for a small number of frequency bins and can be easily implemented in open source software such as Octave or Python.

More info is available on Wikipedia https://en.wikipedia.org/wiki/Goertzel_algorithm including the handy rule of thumb as to when it is more efficient:

$$M \le \frac{5N_2}{6N}log_2(N_2)$$

Estimating A Single Tone in Low SNR Conditions

For estimation of a single tone in low SNR conditions see this paper by Rim Elasmi_Ksibi, Hichem Bessbes, Roberto Lopez-Lacarace and Sofiane Cherife: https://www.researchgate.net/publication/220228033_Frequency_estimation_of_real-valued_single-tone_in_colored_noise_using_multiple_autocorrelation_lags which extracts the estimate for $cos(\omega_0)$ from the samples of the autocorrelation, where $\omega_0 \in [0, \pi] $.  Which derives the estimate for $\omega_0$ as:

$$\cos(\hat{\omega_0}) = \frac{\sum_{k=p}^q \hat{r_k}(\hat{r}_{k-1}+\hat{r}_{k+1})}{2\sum_{k=p}^q \hat{r}_k^2}$$

Where $\hat{r}_k$ is the unbiased autocorrelation for the observed samples y given as:

$$\hat{r}_k = \frac{1}{N-k}\sum_{n=k+1}^N y_n y_{n-k}$$

and p and q are any integers q>p large enough such that the noise in the samples compared are independent.  The larger the range of q-p the more processing required but the lower the noise in the estimate (so you can make that trade). If you choose a p that is less than the lag for noise independence, then you will have no advantage from additional processing for p+n until p+n is at the lag where the noise is independent. You can access this from the autocorrelation of the noise process alone to determine the lag at which the autocorrelation is 0. For example with white noise the autocorrelation = 0 for any offset, meaning all noise samples are indpendent in which case p can be as low as 1.

Intuitive Explanation of Autocorrelation As a Frequency Discriminator

The above gives the actual frequency estimate with all scaling accounted for, and you can trade the computational complexity with the noise of the estimate approaching the CRLB as detailed in the referenced paper. What follows is to provide and intuitive understanding of how this works. This works on the simple principle that the product of a sinusoid and a phase shifted version of the same sinusoid is scaled by the cosine of the phase as given by the following trigonometric identity:

$$cos(\alpha)cos(\beta) = cos(\alpha+\beta) + cos(\alpha-\beta)$$

So when the frequency is the same and only the phase is different, the product is:

$$cos(\omega_c t+\phi)cos(\omega_c t) = cos(\phi) + cos(2\omega t + \phi)$$

If we average (low pass filter) the above, the $cos(2\omega t+ \phi)$ term goes to zero and we are left with $\cos(\phi)$. This shows how the product is a phase detector. When we delay and multiply (as done in the autocorrelation!), the delay produces a signal with the same frequency but a phase shift. The resulting phase measured by the phase detector (product) is the change in phase over that delay which by definition is frequency! (Frequency is the derivative of phase).

A commmon frequency discriminator topology is to delay and multiply (a frequency discriminator produces an output value that is proportional to the frequency of the input):

Each sample of the Autocorrelation Function is a delay and multiply with a different delay value for each. The above referenced paper is simply scaling each result back to be $cos(\omega)$ and averaging to minimize the noise contribution and improve the estimate. In the plot below the vertical axis is crossing the horizontal axis at $-\pi/2$ to be at the point of maximum slope (operating point when used as a discriminator):

This all applies to complex signals as well, in which case a complex conjugate multiply is done as shown with phase detector topologies below comparing real signals to complex signals. This suggests for a single complex tone the use of either the real or imaginary output of the complex conjugate multiplication to get a similar $cos(\omega_0)$ (real out, I) or $sin(\omega_0)$ (imag our, Q) result but with further processing a direct result of $\omega_0$ is obtained using:

$$\omega_0=atan2(Q, I)$$

Where atan2 is the 2-argument arctangent with I and Q are the real and imaginary results of the complex conjugate multiplication, suggesting how the referenced approach for a single real sinusoid can also be extended to the case for a single complex tone.

And for a single complex tone in high SNR conditions the estimate is trivial since the normalized angular frequency is the phase change from one complex sample to next, which is readily extracted from complex conjugate multiplication of the two samples:

$$Ae^{j\omega_0} = y[n-1]y[n]^*$$

With $\omega_0$ extracted using the atan2 function on the real (I) and imaginary (Q) result of the product as $atan2(Q,I)$. This ends up with the following in terms of $y[n]=I[n]+jQ[n]$:

$$\omega_0=atan2 \bigg( \frac{I[n]I[n-1]+Q[n]Q[n-1]}{I[n]Q[n-1]-Q[n]I[n-1]}\bigg)$$

(And there are numerous efficient estimators for the atan2 process that can be used to further simplify this, including the iterative CORDIC rotator when cycle times to iterate are  more available than multipliers and look up tables.)

What is useful and interesting from this is the imaginary portion of the autocorrelation function for any waveform will be proportional to the frequency offset of that waveform, which is useful for carrier recovery implementations for radio receivers! This is demonstrated below in the result for autocorrelation of a complex additive white Gaussian noise signal with a frequency offset in one direction ($e^{j\omega_o t}$) as plotted on a complex plane showing the real and imaginary terms of the autocorrelation.


dealing with Huge Length (Order) FIR filters I'm studying a project.

My signal is audio 8KHz Mono, 8 bits in Low frequencies 70Hz until 1200Hz.

When I perform the calculus with the transition band is very narrow (Transition Frequency) my Order results are (according to proposed window):

For this Example the worst case:

DeltaFrequency = 1.060180596543088 Hz
dNormalizedDeltaFreq= 1.060180596543088/8000

How is calculated the Order? My Answer is below, check the PD!

Rectangular: 6793
Hamming:24903
Hanning:23393
Kaiser:27353
Blackman:41503

If I try to use 16KHz, before orders will be doubled!

Then, I want to use Hardware only to perform the filtering FIR process, and I'm worry about of Time needed to perform this functionalities.

COMMENT1:

Reduce Sampling Frequency is not optional, I need to deal with 1200Hz signals too!

PD:  (This is not question!)

    if (windowType.toLowerCase().equals("rectangular")) {
      Size = (int)Math.ceil(0.9/dNormalizedDeltaFreq);
    }
    if (windowType.toLowerCase().equals("hamming")) {
      Size = (int)Math.ceil(3.3/dNormalizedDeltaFreq);
    }
    if (windowType.toLowerCase().equals("hanning")) {
      Size = (int)Math.ceil(3.1/dNormalizedDeltaFreq);
    }
    if (windowType.toLowerCase().equals("kaiser")) {
//      Mitra's book Handbook for Digital Signal Processing quotes Kaiser with a
//      simple estimate: N = [-20*log10(sqrt(dp*ds)) - 13]/[14.6(ws - wp)/2pi]
       double Att = 60.0;
      Size = (int)Math.ceil((Att - 7.95)/(14.36*dNormalizedDeltaFreq));
    }
    if (windowType.toLowerCase().equals("blackman")) {
      Size = (int)Math.ceil(5.5/dNormalizedDeltaFreq);
    }

QUESTION 1: it is feasible to use a FIR filter of this Order?

QUESTION 2: Which single device (hardware) help me perform this function?

QUESTION 3: How could be calculated the delay time?

There are three different implementations of FIR filters that I am aware of.

The first is the basic "transversal" FIR in which you sum up the terms of the summation directly. It gets very expensive if $L$, the FIR length, gets long.

The second is the so-called "fast convolution" in which the FIR impulse response is zero-padded to a longer length and the DFT is applied to it once. Then there are two slightly different procedures, "overlap-add" (OLA) or "overlap-save" (OLS, sometimes called "overlap-scrap" which I think is more descriptive) to process blocks of the input into blocks of the output.

Both of these methods are well covered in textbooks.

A third method, not really covered in textbooks at all, is the "Truncated IIR" filters (TIIR) and this might apply in your case.  The simplest example of a TIIR filter is the moving-sum or moving-average that you might find in cascade-integrate-comb (CIC) filters. They have recursion in them, but they are not, mathematically, IIR. Their IIR impulse response gets cancelled or truncated.

If you're trying to implement a long FIR that is a sliding rectangular window (that's a moving sum) or a sliding Hann window or sliding Hamming window, you can do that very efficiently with a TIIR.  A sliding Blackman is a little harder and I don't think TIIR will work for a sliding Kaiser.

I have to go now, but I will return with some references. Or send me an email (it's easy to find) and I will send you copies of the references I have including a "how-to" paper I wrote about how to do biquad TIIR filters in general.


Is there a faster algorithm about this type of FFT my problem is:

A real sequence $C$, the length is  $N$ is multiple of 32, the pattern is as below:

$$C =\left(a_0,\, a_1,\, \ldots,\, a_{\frac N2-1},\,\, b_0,\, b_1,\,\ldots,\,  b_{\frac N2 -1}\right)$$

the other two sequence is $A$, $B$, length is also  $N$, as below:

$$\begin{align}
A &=(a_0,& a_1,& \ldots,& a_{\frac N2-1},&&0, &&\ldots&0&)\\
B &=(0,&& \ldots,& 0,& b_0,& b_1,&\ldots,&& b_{\frac N2-1}&)
\end{align}$$

you can see, $C= A+B$.

the DFT of sequence $C$ is known, but sequences $A$, $B$ and $C$ are all unknown.

I just want to know the DFT of sequence $A$ and $B$.

I know the ordinary method is

perform a backward FFT to get sequence $C$, then
sequence $A$ and $B$ becomes known, then
perform forward FFT about $A$, the DFT of $B$ can be get by subtract DFT
of $A$ from DFT of $C$;

but the computaton is too complicated,

Is there a faster computation method?

Alas, a multiplication by a rectangular window in the time domain is a convolution in the frequency domain.  Thus (assuming no other constraints on or knowledge about A or B), the DFT of A (or B) would have to be the convolution of the DFT of C with a twisty Sinc function, which is computationally far more expensive than 2 FFTs (except by fast convolution, which would require the same number of FFTs as your "ordinary" method).


Artificial Neural Network Example Problem I was looking at an example of a NN but I dont understand in the Answer i where the .68 come from or .6637

This is the link of the tutorial

The tutorial states that each neuron uses the sigmoid activation function:

That funciton is:
$$
\sigma(x) = \frac{1}{1+ e^{-x}}
$$
so just substitute $x = 0.755$ into that to get 0.68026719669.


Any experiences for plotting a stationary wavelet transform? I am experimenting with wavelets for my thesis and am currently working with the stationary WT pywavelets provides.
There are very nice plots for CWTs, but does anyone know a technique for producing a plot that gives a good overall understanding of the produced SWT?
Right now I am basically just producing a list of details coefficient plots for each level, and not regarding the averages at all.

I hope this question belongs here.
Cheers!

EDIT:

I really just iterate through the transform and plot each details coefficient vector.
(I'm aware the example uses DWT and not SWT, but it should be analogous)

fig, axes = plt.subplots(3, 5, figsize=[14, 8])

c = pywt.wavedec(data1, 'haar', mode='periodization')

for i in range(0, 5):
        axes[0, i].plot(c[i],c = "r")

The result then looks like this:

I don't do Python, I'm an old person sticking to his old Matlab (codes and) habits. However, up to extension/wavelet/border issues,  SWT is a discrete equivalent to CWT. And in most versions, the number of samples in approximations or details is the same (which is not the case for the DWT). Hence, you can concatenate 1D rows of details (or approximations) into 2D images, akin to traditional scalograms. The following images and code show the process. I have generated a random piecewise polynomial signal with increasing degrees. Then rows of details, and an image of the rows concatenated. The motivation is to address the impact of vanishing moments on piece-wise polynomials, with border effects.
Here are two different realizations.

The Matlab code, that may be reproduced in Python:

dataLength = 512;
data = zeros(dataLength,1);
nChunk = 4;
lChunk = dataLength/nChunk;
for iChunk = 0:nChunk-1
    idxChunk = iChunk*lChunk+(1:lChunk)';
    polyChunk = rand(iChunk+1,1)-0.5;
    dataChunk = polyval(polyChunk,linspace(-0.5,0.5,lChunk));
    data(idxChunk) = dataChunk;
end

nLevel = 4;
[swa,swd] = swt(data,nLevel,'db2');

figure(1);clf
subplot(1,3,1)
plot(data,'x');axis tight
for iPlot = 1:nLevel
    subplot(nLevel,3,3*(iPlot-1)+2)
    plot(swd(iPlot,:));;axis tight
end
subplot(1,3,3)
imagesc(swd)

Other related questions/answers:

STFT and DWT (Wavelets)
Does SWT/ISWT require intermediate approximation coefficients to represent/reconstruct the original signal?
Extract approximation and detail coefficients
Undecimated DWT vs. CWT. In what cases is one preferable over another?
Additive white Gaussian noise and undecimated DWT
Reproducing paper results about a wavelet transformation using python


Frequency Response FIR Filter Bank I want to determine the frequency response of a filter bank
from the input x[n] to the four output channels.

$H_0$ is an lowpass and $H_1$ is a highpass filter.
Both are FIR filters and I have the numerator coefficients of them,
they have an order of 50.

My attempt in matlab would be to use the downsample command on the filter coefficients and multiply them element-wise and pass the result to the freqz command.

So from $x[n]$ to $x_4[n]$ I'd do it the following way:

freqz(downsample(num_h0, 2) .* downsample(num_h0, 2) .* downsample(num_h0, 2));

Is this the right way or am I missing some point here?

Also how do I plot all four frequency responses in one plot,
respectively how do I have to scale the frequency axis?
Is it right to say that the sample rate of $x_4[n]$ is
$\frac{1}{8}$ of the inputs sample rate.

You mix up the convolution theorem: You cannot multiply the impulse responses, but merely their spectra, if you want to go this way.

However, a simpler way (esp. when you have Matlab available and look for a numeric solution) is to calculate the overall impulse response. For example, considering the filter $x[n] \rightarrow x_1[n]$ you have:

Let $x[n]=\delta[n]$. Then, $x_1[n]=h_1[2n]$, i.e. the impulse response of the overall filter is $h_1[2n]$. Note, that n is zero-based here (i.e. first sample is at $n=0$).

You can also do this in Matlab. Generate a delta impulse and send it through the filter bank. The resulting impulse responses are your effective FIR coefficients for each filter path:

N = 50;
h0 = ones(N, 1).*exp(-(1:N)');            % a lowpass
h1 = sin(3*(1:N)).*exp(-0.3*(1:N));     % a highpass

x = [1; zeros(1,1)];  % a delta impulse

x1 = downsample(conv(x, h1), 2);

x2p = downsample(conv(x, h0), 2);
x2 = downsample(conv(x2p, h1), 2);

x3p = downsample(conv(x2p, h0), 2);
x3 = downsample(conv(x3p, h1), 2);

x4p = downsample(conv(x3p, h0), 2);
x4 = x4p;

% x4 contains the impulse response of the filter x->x4
f = 1/pi*linspace(-pi, pi, 1000);
f1 = 1/pi*linspace(-pi/2, pi/2, 1000);
f2 = 1/pi*linspace(-pi/4, pi/4, 1000);
f3 = 1/pi*linspace(-pi/8, pi/8, 1000);
f4 = 1/pi*linspace(-pi/8, pi/8, 1000);

subplot(2,1,1);
plot(f, mag2db(fftshift(abs(fft([h0, h1'], 1000, 1)), 1)));
xlim([0, 1]);
xlabel('f/pi');

subplot(2,1,2);
hold off;
plot(f1, mag2db(fftshift(abs(fft(x1, 1000)))));
hold on;
plot(f2, mag2db(fftshift(abs(fft(x2, 1000)))), 'r');
plot(f3, mag2db(fftshift(abs(fft(x3, 1000)))), 'k');
plot(f4, mag2db(fftshift(abs(fft(x4, 1000)))), 'g');
hold off;
grid on;
xlabel('f/pi');
xlim([0, 1]);

As you see that the sample rate is reduced for each path, the maximum frequency is also reduced by half. If you were to plot all frequency responses over the full frequency range, you'd need to periodically repeat them.


Lagrange Multipliers Optimization - Complex Functions I want to optimize a sum of complex frequency responses (complex sum is SPL, sound pressure level as a complex function of distance) of a number of loudspeakers, say 10. There is a reference, desired SPL curve given (objective function), which represents the magnitude of the sum of 10 frequency responses. I have read a number of papers how to do optimization using Lagrange Multipliers, however, did not find many papers (only one) about optimization when the functions are complex (Re and Im parts).

Any help or pointing to a good references, examples would be appreciated.

I found the following in Charles Therrien's "Discrete Random Signals and Statistical Signal Processing" in one of the Appendicies.

Say you have the function $Q(a)$ you wish to minimize such that $C(a)=0$, where $C(a)$ may be complex valued and $a$ may be a complex vector. The constraint really represents two real-valued constraints.
$$C_r(a)=0,\qquad
C_i(a)=0,
$$
where $C_r(a)$ and $C_i(a)$ represent the real and imaginary components. The Lagrangian can be written as:
$$\mathcal{L}=Q(a) +\mu_1C_r(a)+\mu_2C_i(a),
$$
where $\mu_1$ and $\mu_2$ are real valued.

This can also be written more compactly using a single complex parameter $\lambda=\lambda_r+j\lambda_i$ as
$$\mathcal{L}=Q(a) +\lambda C(a)+\lambda ^*C^*(a)
$$
where $\lambda^*$  and $C^*(a)$ denote the complex conjugate of the quantities. The reason for this is that
$$
\lambda C(a)+\lambda ^*C^*(a)=2Re[\lambda C(a)]=2\lambda_rC_r(a)-2\lambda_iC_i(a)
$$
and letting $\mu_1=2\lambda_r$ and $\mu_2=-2\lambda_i$, then the two Lagrangians  are seen to be equivalent.

Note - to apply the gradient to a real valued function with respect to a complex parameter you should also be familiar with Brandwood "A complex gradient operator and its application in adaptive array theory" IEE Proceedings 130(1) 11-16, Feb, 1983. This essentially treats the variables $a$ and $a^*$ as two separate variables and you only have to take the derivative wrt to one of them to perform a minimization.


Removing specific frequency from a signal How do you remove a specific frequency from a signal?

I have a signal which is basically a sine wave at 20hz with some randomness added to it. After the real fft I can clearly identify the peak in the 20th bin, but when zeroing this bin and those sorrounding it and converting back to time domain the whole signal is gone for good, where I would expect to see the "random" part of the signal more or less intact.

Below is a picture that has:

Top graph has the noiseless 20Hz sine wave sampled at 100Hz (blue) and the noisy sine wave (red).
Bottom graph has the actual noise signal added (red) and the noisy signal filtered with an IIR notch filter with a notch centred at 20Hz (blue).

Apart from the group delay between the original noise and the output of the notch filter, the two are very similar.

R code to implement this is below the graph.

The frequency response of the notch filter is as below.

R Code

#26866

f <- 20
fs <- 100
phi <- 2*pi*0.12987892374

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

signal <- sin(2*pi*f/fs*t+phi)
noise <-  rnorm(T,0,1)
noisy_signal <- signal + noise

alpha <- 0.9

num <- c(1, -2*cos(2*pi*f/fs), 1)
den <- c(1, -2*alpha*cos(2*pi*f/fs), alpha*alpha)

filtered <- signal::filter(num,den,noisy_signal)

par(mfrow=c(2,1))
plot(noisy_signal, type="l", col="red",  lwd=10)
lines(signal, type="l", col="blue", lwd = 3)
title("Noiseless signal (blue) and noisy signal (red)")

plot(filtered, col="blue", lwd = 4)
lines(noise, col="red",  lwd=3)
title("Filtered signal (blue) and original noise signal (red)")


What is this strange pulse in my audio output? Firstly, if this is the wrong place for this, I'd be happy to know the correct place for this and I'll move it there.
I have two PCs which share an audio interface. My laptop is connected via USB and uses it as its sound card. My Desktop (runs Linux so will be called Linux PC from here on) inputs into inputs 3 and 4 of the interface which is routed out the monitoring outputs. If it's relevant, I'm using a Focusrite Scarlett 6i6 (2nd Gen). I was watching something on YouTube on the Linux PC about an experiment involving sound, one of those things that let you try at home. Interested, I opened my DAW on my laptop, set it to read inputs 3 and 4 (the Linux PC), and opened a spectrum analyser to take a look at the frequencies of sound.
My spectrum analyser shows the audio has a pulsing every 10 seconds with the pulse lasting around 7 seconds, with two pulses in unison at ~20kHz reaching up to ~24kHz and the second at ~40kHz reaching up above the scale which tops off at ~43kHz. My sample rate is 88.2kHz.
I've included a screenshot. The left is a spectrogram before applying a High Pass filter (via EQ) set at around 18kHz, the right is after filtering.

What is this?
I've looked at this question and accepted answer, which is the only thing I could find that has come close to be even remotely similar to what I'm seeing and I don't quite know if it's related or not. The post mentions to filter out anything above 1/2 the sample rate to filter out aliasing artefacts, but 1/2 my sample rate is 44.1kHz which plants the lower pulse well below this area and the higher one just below it.
The screenshot was taken when (seemingly) no audio was playing. I can't tell if the plugin has some kind of automatic gain control that makes the maximum amplitude fit in the spectrum's colour gradient or not, but when audio plays it is significantly weaker on the spectrum but definitely still there.
I've also recorded, intensely amplified and lowered the pitch and it sounds like data, but again, what data?

Not really an answer but some hints on how to proceed.

You need to debug this on step at a time. Start with shortening the inputs to the Focusrite and disconnect the Linux computer entirely. If you still see noise, disconnect the power supply from you laptop and run battery only. If it's still noisy, either the Focusrite or the USB interface are bad (i.e. ground  provided by your USB port is just too dirty).
Then connect the Linux PC ground only (leave the inputs shorted).
Then connect signal as well but make sure the output of the Linux PC is all the way down and verify by measuring: should be zero volts.
Then play a -20dBdFS sine wave at 1kHz from the Linux PC. Compare that to expected results.

The fix will depend on where in the chain things go wrong. Audio connections are tricky and susceptible to power supply contamination, ground loops, Electro Magnetic Interference, etc.
Hence you should used balanced connections whenever possible. Try using the XLR inputs on the Focusrite by either using a properly configured "mono to balanced" adaptor cable. Means that you should NOT connect signal ground. You can play around with connecting chassis ground or not. The easiest solution would probably be to use a passive audio isolation transformer with a ground lift switch.  Something like this perhaps: https://www.audiopile.net/FDB-202 (I'm not affiliated).
Good reads about this: https://www.ranecommercial.com/kb_article.php?article=2107


How does sampling rate affect discrete filters? I am just starting to learn discrete filters and I could use some help. I understand continuous signals and filters.

I am trying to understand the math behind discrete filtering. For example in the s domain, a simple low pass filter can be recognized as $$\dfrac{1}{\tau s+1}$$
 How does sampling rate affect a discrete low pass filter?

IF we take the sampling period as T.First I will convert the filter to time domain impulse response as
$$
h(s)=\frac{1}{\tau s+1}\\h(s)=\frac{1/\tau}{s+\frac{1}{\tau}}
$$
this in time domain is
$$
h(t)=\frac{e^{-t/\tau}}{\tau}
$$
IF the above signal is sampled at time T the it can be written as
$$
h(Tn)=\frac{e^{\frac{-Tn}{\tau}}}{\tau}\\h(n)=\frac{e^{\frac{-Tn}{\tau}}}{\tau}
$$
if Z transform taken for the impulse response.Then the equation is
$$
H(z)=\frac{1/\tau}{1-e^{\frac{T}{\tau}}z^{-1}}
$$
Here we can see the sampling time coming in denominator which could put the pole inside or outside the unit circle


Why can Quadrature Demodulation demod a Frequency Modulated Signal? I use a Quadrature Demodulator in my SDR application, which is defined as:

$\angle (S_n, S_{n+1})=arctan(S_{n+1}*\overline{S_{n}})$

So practically its amplitude is the angle between two Samples $S_n$ and $S_{n+1}$ where $n$ refers to a position in a sequence of complex values of I and Q. $\overline{S}$ be the complex conjugate.

I understand that for Quadrature based modulation the phase-changes contain the information. But why can I demodulate frequency-modulated Signals with a Quadrature Demodulator. This modulator is often referred to as FM Demod. Can somebody explain me why that is?

Best,
Marius

I don't think your formula is quite right.

The two complex samples are the locations (at two successive sampling instants)
of the tip of the rotating phasor that represents the analog signal.  The information that you need is the angle between phasor positions at these
two successive time instants.  If $S_{n+1} = r_{n+1}e^{j\theta_{n+1}}$ and
$S_n = r_ne^{j\theta_n}$, then you want the value of $\theta_{n+1}-\theta_n$.
Now, $S_{n+1}S_n^* = r_{n+1}r_ne^{j(\theta_{n+1}-\theta_n)}$ and if you
express these complex numbers in rectangular coordinates, then you can
get the desired angle as
$$\theta_{n+1}-\theta_n
= \arctan\left(\frac{\text{Im}(S_{n+1}S_n^*)}{\text{Re}(S_{n+1}S_n^*)}\right)
= \arctan\left(\frac{\text{Re}(S_{n})\text{Im}(S_{n+1})
- \text{Im}(S_n)\text{Re}(S_{n+1})}{\text{Re}(S_{n})\text{Re}(S_{n+1})
+ \text{Im}(S_{n})\text{Im}(S_{n+1})}\right)$$
though usually the four-quadrant arctangent function atan2(y,x) is used
instead of atan which returns an answer between $-\pi/2$ and $\pi/2$.

In digital communications applications such as DQPSK demodulation,
the real interest is not in the actual value of the angle, but in which
of the four intervals $(-\pi/4,\pi/4)$, $(\pi/4, 3\pi/4)$, $(3\pi/4, 5\pi/4)$,
and $(5\pi/4, 7\pi/4)$ the angle belongs -- in other words,
did the phase not change at all, or change by $\pi/2$, or by $\pi$, or
by $3\pi/2$ -- and this can be readily determined
by comparing the values of the numerator and denominator in
the argument of the arctangent above, thus saving computational
effort and speeding up the demodulation process. But the
more general formula gives the actual angle, and thus can be
used to demodulate a general FM signal as well.  The output
will be a sequence of samples of the original continuous-time
signal that modulated the carrier and thus formed the signal
that was transmitted.  D/A conversion will give the reconstructed
signal.


Using ROC curves for comparing the performance of SIFT and SURF I want to compare the performance of SIFT and SURF against some standard image set.

After a bit of research, I found that ROC curves are one of the best ways to do this task.

In order to plot ROC curve TP rate and FP rate are needed.

Output of SIFT and SURF are detected keypoints. How exactly can I generate ROC curves from those keypoints ?

The ROC:
ROC curves are popularly used as performance metric for classification tasks. If the images in your dataset has class labels, then you can employ supervised learning to train a classifier (SVM for example). The dataset is split into training and testing and predicted class score from the classifier for images in the test set are compared to ground truth labels. The ROC curve is acquired by applying a threshold value to the classifier predicted score and obtaining a (TP,FP) value for each threshold to generate the curve.

The features:
SIFT and SURF are detectors and image local-patch descriptors (you can also use these without the detector in what is now more successful approach called 'dense sampling' where the keypoints are uniformly sampled in the image).

There are several approaches, but the most popular and recommended for first use is the 'Bag-of-Words' (BoW) model based approach. The reason is that: each image has one label (ignoring multi-class and multi-label scenario for now) and several hundred keypoints and each keypoint has an associated high-dimensional feature vector. The feature for an image at this stage is a matrix and we would like a vector (a point in feature space for each image with its label). The BoW approach does just this job for us.

BoW approach:
All the feature vectors for every image in the training set are collated and clustered (ignoring labels and location data) to compute a 'dictionary', using unsupervised learning. The dictionary is a set of cluster centroids, most commonly computed using k-means clustering. The dictionary follows data compression scheme to succinctly represent the data distribution.

Subsequently, each image is 'encoded' using its keypoints and the computed dictionary. Using vector quantization VQ, each keypoint is associated with a dictionary element to compute a occurrence histogram of keypoint associations. The consequence is that each image now has a feature vector which is the normalized occurrence histogram of keypoints assigned to dictionary elements. The dimensionality of this encoded feature vector is the size of the dictionary.

This normalized histogram for each image along with the label is used to training a classifier and then acquire classification results measured in terms of ROC curve. The sequence is:

[Images] -> (SIFT/SURF) -> [SIFT/SURF feature matrix] -> (collate and cluster) -> [dictionary] -> (BoW histogram by VQ encoding) -> [feature vector, label] -> Classifier training (training data) -> Classifier testing (validation/test data) -> ROC curve.

Suggestions:
Besides ROC, there are other performance measures available like accuracy, mean average precision, F1 score, etc. Your choice depends on the imageset, label ratio balance, objective of the classifier, etc.

There is plenty of code available online for baseline Bag-of-words model for images. One sample in Matlab using the Caltech101 dataset is at:
http://www.vlfeat.org/applications/apps.html#apps.caltech-101

The performance of SIFT and SURF is also provided in the original paper on SURF, so you can compare your results if you choose to work with the same dataset.


Low data rate with high frequency offset I am trying to implement a basic digital modulation scheme (BPSK) with a low data rate of 5 kbps. The requirement is to support a frequency offset of +/- 10 KHz.
I would like to know what are the criteria to decide how much offset can be corrected? My burst data header is only 50 bits. I want to know if it is possible to track the frequency offset and also do the timing recovery (gardner) in the header bits.

The answer is yes it is possible and any amount of frequency offset can be corrected to the extent your sampling rate is sufficiently high enough for unambiguous frequency resolution over ±10 KHz. A very simple frequency discriminator that can be used in a carrier tracking loop is further described at this post:Carrier frequency offset estimation for QPSK and asymmetrical spectrum and results in determining the frequency offset simply from the product of two successive (complex) samples given as $I_1Q_2-I_2Q_1$. For BPSK you can also remove the modulation from the frequency error result by first squaring the signal, in which case the result will be twice the actual frequency offset; which is of no consequence since the loop will drive the offset to 0 in any event).
This discriminator has unambiguous frequency resolution for IQ rotation up to 90° between samples. A 10 KHz frequency offset will rotate 90° if the sampling rate is 40 KHz, so as long as your sampling rate can be greater than this rate, this should be a very simple approach to carrier recovery / Doppler wipe-off (80 KHz if using the frequency doubling approach to strip the modulation).
The maximum possible offset sets your search window and the header duration if used for acquisition sets the carrier recovery Loop BW. The carrier Loop BW sets the acquisition SNR required (wider BW is faster but requires a higher SNR). So use a loop BW that is wide enough to settle for a step response in frequency change from 0 to 10 KHz within the time duration of 50 symbols. (And then the loop bandwidth can be further tightened once acquired if channel dynamics allow).


Direct form 1 filter equation for a biquad filter In general, we all know that the direct form 1 equation for a biquad filter is given by
$$y(n) = b_0x(n) + b_1x(n-1) + b_2x(n-2) - a_1y(n-1) - a_2y(n-2)$$
But when it comes to practical applications like for audio processing applications, I have seen many people using
$$y(n) = b_0x(n) + 2b_1x(n-1) + b_2x(n-2) - 2a_1y(n-1) - a_2y(n-2)$$
as the filter equation in their codes.
So, I just want to know the reason why we need to multiply $a_1$ and $b_1$ coefficients by 2?

This may be a misunderstanding based on a fixed point implementation. Fixed point implementations can benefit from "fractional multiplies" which requires a number format where all absolute values are  1 or less. In Q notation that would be Q15.0 for a signed 16-bit integer. Many processors have a dedicated hardware instruction for this operation, which makes it fast and efficient.
Assuming that all poles and zeros are inside the unit circle, we have $|b_0|, |b_2|, |a_2| < 1$ but $|b_1|, |a_1| < 2$. In order to avoid shifting associated with mixed Q notation, the $b_1$ and $a_1$ coefficients are split in have, so numbers can be kept in the same Q. Hence it would be implemented as
$$y[n] = b_0 \cdot x[n] + b_1/2 \cdot x[n-1] + b_1/2 \cdot x[n-1] + b_2 \cdot x[n-1] ...$$
or
$$y[n] = b_0 \cdot x[n] + (b_1/2 \cdot x[n-1]) \ll 1 + b_2 \cdot x[n-1] ...$$
where $\ll$ is the left-shift operator.
This may looks like the coefficient is multiplied by two but in reality it's first split in half to make it's magnitude smaller than $1$, so the same type of fixed point multiplication can be used for all operations.


Artifacts signal conversion from float to integer I have a script where i want to change the noise floor to about 96 dB (visual purpose).
The script containing basically a sine wave with after wart do FFT over it.

The script is in python and as follows:

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import get_window
from scipy.fftpack import fft

l = 128
f = 4
x = np.arange(l) / l
max_int = np.iinfo(np.int16).max
sine = np.sin(2 * np.pi * f * x)
sine_96 = (sine / sine.max() * max_int).astype(np.int16)

w = get_window('hann', l)
S = np.abs(fft(sine * w))
S_96 = np.abs(fft(sine_96 * w))

fig, ax = plt.subplots(2, 2, figsize=(8, 7))
ax[0, 0].plot(sine * w)
ax[0, 1].plot(sine_96 * w)

ax[1, 0].plot(10 * np.log10(S[:l // 2]))
ax[1, 1].plot(10 * np.log10(S_96[:l // 2]))

The weird thing is when i change the signal to 16 bit integer (or 16 bit integer and back to floats again) i get some artefacts what i don't understand. See also

left are the floating point functions where i get about an 150 dB noise floor. and on the right the version where the signal is changed to 16 bit integer. The result should be the same signal but with a noise floor around 96 dB.

As i can see it the noise floor is getting higher but so do the peak. I also get's repetitive peaks.
I have two question how can i avoid these peaks and what is the background of this?

The quantization error is dependent on the input signal. In your case the input signal is a sine wave and the quantization error is a combination of its 8 odd harmonics that fit into the available bandwidth. You can reduce the dependency of the quantization error to the input signal by adding triangular dither before quantizing: Add to each sample a uniformly distributed random number between 0 and the quantization step, and subtract from the result another uniformly distributed random number between 0 and the quantization step. Then quantize. This way you'll get a flat spectral noise floor (no distortion of the signal) and no modulation of the quantization noise by the signal.

For further reading on dither: Robert Alexander Wannamaker's PhD thesis "The Theory of Dithered Quantization", 2003.

Modifying your script to include dithering fixes the problem:

import numpy as np
import random
import matplotlib.pyplot as plt
from scipy.signal import get_window
from scipy.fftpack import fft

l = 128
f = 4
x = np.arange(l) / float(l)
max_int = np.iinfo(np.int16).max
sine = np.sin(2 * np.pi * f * x)
dither = [random.triangular(-1, 1, 0) for _ in xrange(l)]
sine_96 = (sine / sine.max() * (max_int - 1) + dither).astype(np.int16)

w = get_window('hann', l)
S = np.abs(fft(sine * w))
S_96 = np.abs(fft(sine_96 * w))

fig, ax = plt.subplots(2, 2, figsize=(8, 7))
ax[0, 0].plot(sine * w)
ax[0, 1].plot(sine_96 * w)

ax[1, 0].plot(20 * np.log10(S[:l // 2]))
ax[1, 1].plot(20 * np.log10(S_96[:l // 2]))
fig.show()

Figure 1. Top left: windowed input signal, top right: dithered, quantized, and windowed input signal, bottom left: frequency spectrum of the windowed input signal (vertical axis in dB scale), bottom right: frequency spectrum of the dithered, quantized, and windowed input signal (vertical axis in dB scale).

Let's plot the quantization residual (error) of the undithered case before windowing, and the frequency spectrum of the residual after windowing it:

sine_96 = (sine / sine.max() * max_int).astype(np.int16)
residual = sine_96 - (sine / sine.max() * max_int)
fig, ax = plt.subplots(1, 2, figsize=(8, 3))
ax[0].plot(residual)
R = np.abs(fft(residual * w))
ax[1].plot(20 * np.log10(R[:l // 2]))
fig.show()

Figure 2. Quantization residual and its spectrum (window function applied before transforming).

Quantization is equivalent to a memoryless waveshaper shaped like a staircase. For a periodic input signal a memoryless waveshaper can only produce harmonics of the input signal, because the resulting output signal will be identically periodic:

$$f(x + s) = f(x) \text{ for all } x \in \mathbf{R}\\
\Rightarrow g\big(f(x + s)\big) = g\big(f(x)\big) \text{ for all } x \in \mathbf{R}$$

Here the input signal is a sine wave with a period of 128/4 = 32 samples. If the period would not be an integer number of samples, then also other frequencies than the harmonics would be generated, because quantization operates in discrete time and non-integer periodicities don't exist in discrete time.


Extract a signal that is known to repeat every 4 samples at a known phase i.e. exactly 500 Hz signal in 2000 Hz data I have a situation that the measurement device is generating a known signal in the data at 1/4 of the rate. In this case that is a 500 Hz signal when sampling at 2000 Hz.

Now I can extract the 500 Hz signal using a Butterworth Band-Pass at 490–510 Hz which works, but it takes >1/8th of a second to extract the signal. In our setup we know the measurement and the signal are perfect in sync (same clock source etc) so a perfect 500 Hz filter would be ideal.

Is there a better extraction method? Ideally without the latency of the Butterworth or more computationally efficient or more accurate to extract the signal?

Here is an example of the 2 kHz data:

-34561, -32650, -3422, -5064, -33967, -32061, -2807, -4436, -33364, -31447, -2108, -3771, -32863, -31050, -1801, -3516, -32634, -30871, -1711, -3496, -32761, -31115, -2027, -3919, -33237, -31652, -2684, -4606, -33831, -32279, -3394, -5291, -34442, -32775, -3764, -5582, -34667, -32969, -3872, -5589, -34532, -32665, -3446, -5064, -34014, -32082, -2771, -4385, -33355, -31437, -2114, -3771, -32854, -31032, -1796, -3501, -32651, -30819, -1593, -3425, -32671, -31067, -1924, -3859, -33214, -31631, -2639, -4519, -33802, -32225, -3301, -5254, -34406, -32736, -3741, -5594, -34652, -32932, -3886, -5572, -34453, -32592, -3377, -5056, -33953, -32043

Edit:
 The aspect of interest is the amplitude of the 500hz signal, the amplitude can be converted to give a circuit impedance which we wish to efficiently measure. The butterworth takes a long time for a 1st order butterworth to converge to a steady amplitude measurement with fixed circuit impedance. A higher order could viable but this feels like overkill and as there is so many known about the signal so a more efficient approach would feel to be viable, just magnitude is the key. I tried some simple approached but these still had oscillation of amplitude at 50hz albeit the average was good, there may be formal methods I'd like to try.

Update: Here is a plot showing the butterworth vs the given solution shown in Green which does exactly what I had hoped with 4 samples latency! Previously I had looked at peak-signal but taking peak-average does look preferable anyway and what the solution is giving us. Thankyou for your help.

To estimate the amplitude of the 500 Hz signal, a simple moving average filter will give you the best estimate in the presence of white Gaussian noise. This can be easily done by doing a moving average on the result of converting the 500 Hz signal to baseband as described below. This filter will have a Sinc response in frequency, so has high sidelobes which means it would be sensitive to the presence of other stronger signals if present. A very simple work-around is to just window the time domain signal first and then follow the same process (downcovert to baseband and sum  the result over N samples where N is the length of the moving average).   This is essentially a bandpass filter (using the windowing filter design approach) but in any filter used there will be a trade between frequency selectivity and time of convergence.

A convenience of having a "digital IF" that is exactly $Fs/4$ is you can down-convert to baseband (create the equivalent baseband analytic signal) by just multiplying by +/-1 as follows:

To get I you multiply by

I = 1 0 -1 0 1 0 -1 0

To get Q you multiply by

Q = 0 -1 0 1 0 -1 0 1

As these would be the values for the Local Oscillator $e^{-j\omega t}$ for $\omega = \pi/2$ which is the normalized radian frequency at an LO of $f_s/4$.

I detail this approach further at this post:

DAC and ADC architecture in SDRs


Prewitt operator and central difference? When taking the gradient of an image one can do $I * \begin{bmatrix}+1 \ 0 \ -1\end{bmatrix}$ (x direction).

As I understand this is basically applying central differences in the x direction, $\delta_h[f](x) = f(x+\tfrac12h)-f(x-\tfrac12h)$. But since the gradient is a vector of the derivatives and $f'(x) \approx \frac{\delta_h[f](x)}{h}$ shouldn't it be $I * \frac{1}{2} \begin{bmatrix}+1 \ 0 \ -1\end{bmatrix}$?

And what about the the Prewitt operator $\begin{bmatrix}+1 \ 0 \ -1 \\+ 1 \ 0 \ -1 \\ +1 \ 0 \ -1\end{bmatrix}$ (x direction)? As I understand there's some relation to finite differences as well, but how does having a 3x3 kernel affect this relation?

Yes, the $1/2$ factor correction could be present, so that the magnitudes between a) the continuous derivative and b) the approximated gradient remain consistent in some way: a (continuous) line with a unit slope will have its discretized version get a unit gradient with the $1/2$ factor.

However, as long as the factor is nowhere applied, all gradient magnitudes remain comparable in a relative way. I suspect that, Prewitt being an ancient operator, computations were of concern. As is only involves simple additions and subtractions, one probably did not care about "costly" divisions, and most of all returning non-integer values.

The same happens with other kernels, see using a sobel operator - edge detection. The normalization is somewhat arbitrary, using a sum of absolute values, or in energy, as for Fourier transform normalizations.

However, if one uses different kernels, or different orders of derivations, it is safer to perform a consistent normalization, in order to be able to compare magnitudes.

The idea behind the constant columns is that a vertical edge should spread over several pixels, so you perform an additional averaging in the direction parallel to the edge, to reduce noise effects for instance, providing you with a more robust estimate of the actual gradient.


Frequency change detection in step-wise frequency chirp I am driving a system with a step-wise frequency chirp from 50-400 hz with a frequency step of 0.01 hz. I am having difficulty identifying the points of frequency change as the data collected has an indeterminate amount of time at 50 Hz and an indeterminate amount of time at 400 Hz. Since the time between frequency changes is 0.1 seconds an fft calculation doesn't provide the proper frequency resolution to determine this change. I also tried calculating the gradient of the frequency chirp with inconclusive results as well as attempting a lock-in calculation. I wrote some sample code that has a randomly located frequency step followed by another 0.1 secs later that I have been trying to diagnose with varying success. Ultimately I am trying to identify the beginning and end of the frequency chirp since at that point the step changes are known but this change point identification was attempted to identify the point of change between 50 and 50.01, then backtracking to find the beginning of the series. Any suggestions would be very much appreciated.

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

fs = 10000
t = np.linspace(0,1,1*fs)
x = np.zeros(len(t))

step_loc = int(2500 + 5000*np.random.random(1))

x[:step_loc] = np.sin(50*2*pi*t[:step_loc])
x[step_loc:step_loc+1000] = np.sin(50.01*2*pi*t[step_loc:step_loc+1000])
x[step_loc+1000:] = np.sin(50.02*2*pi*t[step_loc+1000:])

grad = np.gradient(x,t[1])

lock = x*(np.cos(50*2*pi*t) + 1j*np.sin(50*2*pi*t))
b, a = signal.butter(5,30/10000, 'low')
lock_filt = signal.filtfilt(b, a, lock)

plt.plot(t,x)
plt.plot(t,grad*3e-3)     #this is offering very inconsistent identifiable results
plt.plot(t,lock_filt*1e1) #although this can identify the general area of the random frequency change i
                      # dont know how to extend it to identifying the actual point

I suggest monitoring the phase versus time directly instead of frequency. Frequency is the derivative of phase so the slope of the phase would indicate the frequency. Detrend the phase slope for the starting frequency and then the point in time where the phase starts to ramp up should be easier to detect. The window in which to detect this change will be balanced with the SNR of the signal itself.

Since frequency is the derivative of phase, it will be more sensitive to high frequency noise (the derivative is a high pass function), which is why monitoring phase directly could be a more reliable approach.


Unstable frequency/phase when using external reference for Ettus B200 I am attempting to do a simple record of a CW, 250MHz signal at the Rx port of an Ettus B200 using an Anritsu MG369 signal generator. I am using the 10 MHz ref signal provided by the SG (output @ ~9dBm). The frequency/phase of the signal seems to be too unstable from what I should expect. Is this common for the B200? I am expecting a B210 soon and am just trying to get ahead of the game while awaiting the new hardware, but if this is user error I'd like to figure it out in advance if possible as I eventually will need to consider phase closely when test is in its final form. I've seen in Ettus literature that a square wave is best for reference, could a sine reference be the cause of this?
Could it be the PLL misbehaving or an amp somewhere in ref circuitry? RF front end is confirmed to be working okay and can see a clean, stable signal when running CW from Tx through cable to Rx port.
Including flow chart as well as two pics of spectrum to try to illustrate instability.
EDIT: Including better FFT sink snapshots

You're looking only at the real part. Obviously, if the phase of transmitter and receiver aren't locked (like here), that will mean the sent oscillation will be more or less in-phase.


Where in the signal chain is Channel Equalization applied? I am reading about channel equalization and I am somewhat confused on where in the signal chain channel equalization is normally applied and if it's before or after a correlator detector or matched filter.

My confusion is that if equalization occurs before matched filtering/correlator detection, then a training sequence used to estimate the channel will need to be based upon the output of the ADC, not of the matched filter. So if a pseudo-random sequence is used as a pilot tone it won't consist of symbols but of some other modulation type.

If the equalization is based upon the symbol rate output of the matched filter how are fractional delays handled?

Generally, any matched filtering would be done before equalization, since there is no practical advantage that I am aware of to applying a matched filter afterwards.  Most textbooks that I've read present an architecture in which a matched filter is applied first, the signal is sampled at 1 sample/symbol, then an equalizer is applied at 1 sample/symbol.  But it does not have to be this way: you could apply a matched filter, but not downsample, then apply a fractionally-spaced equalizer and sample the output at 1 sample/symbol; or your could skip the matched filter entirely and apply an equalizer directly, in which case the equalizer would adapt to include both the matched filter and equalizer response.

Regarding fractional delays: if the equalizer is symbol-spaced, then we must perform symbol timing recovery prior to downsampling so that we sample the matched filter output at the correct times.  We need to perform symbol timing recovery anyway to correct any rate offsets, but if we are downsampling to 1 sample/symbol, we must also get the symbol clock phase correct.  If the equalizer is fractionally-spaced, then it can adaptively refine the sample time, although it can't correct any substantial rate offsets.  This is an advantage of using a fractionally-spaced equalizer; the primary disadvantage is the higher input sampling rate.


area under the phase encoding gradient in MRI I have a question about how I could compute the area (duration and amplitude) for the phase encoding gradient in a typical MRI gradient echo sequence. So, let us assume that my field of view in the phase-encode direction is 20 cms and I have a k-space matrix of size 256. From the tutorial here, I can see that the k-space spacing is $1/20e-3$ per meter. Let us also assume that the maximum gradient strength I can achieve is 10mT/m.

I am now stuck and not sure how I can compute the area that I need under the phase encode area. I am aware that the area will change at each phase encoding step and at least from the pulse sequence diagrams that I see for a GE sequence, it seems that the amplitude is adjusted (going for negative max to positive max) and the gradient time is kept constant during each phase encode step.

However, I am at a loss as to how to compute this area?

Actually, it is as simple as the frequency encoding gradient. Both encoding gradients look terribly different, but in the end they work the very same way.

The outer $k$-space position $k_\textrm{max}$ is related to your resolution, e.g. $1\,\textrm{mm}$, or in terms of the FOV: $FOV/N$ where $N$ is the number of image pixels.
For the steps between two $k$-space lines, $\Delta k = 1/FOV$

With $k(t) = \frac{\gamma}{2\pi}\int_0^t G(t^\prime)\,\textrm{d}t^\prime$ you can find that for the largest gradient moment $M_\textrm{max} = \frac{\pi(N-1)}{\gamma FOV}$. The steps in between scale this down linearly.

Pay attention to the line numbering, because the central line in $k$-space should be acquired with zero gradient moment.


Magnitude characteristic of H(s) I have a question about drawing a magnitude characteristics of H(s).

So my transmittance function is for example:

and in the book the magnitude characteristics is shown as:

And my question is - how to count the value of abs(H(jΏ) ?

I understand that s = 1 is a zeros (local minima), and this function has two poles: -1-j and -1 +1 (local maxima), but I can't count the values highlighted in yellow.

The most direct way is to turn your transfer function into frequency response by setting $$s = j\Omega$$ and then $$H(j\Omega) = \frac{j\Omega - 1}{-\Omega^2 + 2j\Omega + 2}$$ So for $\Omega = 0$, you get $$H(0) = -\frac{1}{2}$$ which gives $|H(0)| = \frac{1}{2}$ and is in agreement with your figure. Similarly, for $\Omega = \pm 1$, you get $$H(\pm j1) = \frac{\pm j - 1}{-1 \pm 2j + 2} = \frac{\pm j -1}{1 \pm 2j}$$ You can show that $|H(\pm j1)| = \sqrt{2/5}$, which also agrees with your magnitude response plot.

However, these are only some interesting values of the magnitude response. If you need to draw the magnitude response in greater detail without using some software to do this for you, you have to go for the geometric approach. Here is a good place to start.


64 QAM versus lower order modulation This question is regarding digital modulation.

If 64 QAM always provides rates greater than BPSK or 16 QAM, then why do systems always use higher order modulation of QAM?

Thanks

There is a trade-off between modulation order, required energy and bit-error rate (BER). As you increase the modulation order, you need to increase the average energy per bit to keep the same BER.

As an easy example, consider BPSK vs 4-PAM. Say you have BPSK with distance between constellation points $d_{min}=1$; this works out to an average bit energy of 0.25 joules per bit.

Now, you want to double the bit rate by going to 4-PAM, but you want to keep roughly the same BER, which means having the same minimum distance $d_{min}$ between constellation neighbors. Now, your average symbol energy is $$\overline{E_s}=\frac{2\cdot1.5^2+2\cdot0.5^2}{4}=1.25$$ joules, or an average bit energy $\overline{E_b}=0.625$.

As you increase the modulation order, you'll see the energy required increase exponentially. Exactly the same thing happens when using QAM.


Hw does a digital piano work? I know someohow it combines different harmonics to synthsis note, but I don't exactly know how this happens. Does digital piano use pwm to synthesis a note?

Piano sound generation is complicated by the sound being so familiar to most of us, and the slowly decaying non-stationary non-harmonic nature of overtones meaning that simple attack + looped samples is hard to get right. It is simplified by the low dimensionality of user interaction (88 keys with velocity, 3 pedals) compared to, say, a guitar where there are many dimensions to how the string is picked and how your left hand interacts with the string.
It used to be that digital pianos used synthesis in order to produce tones. My old Roland piano used «SA synthesis»:
http://www.novelmusic.com/roland/Technical/default.htm
As digital memory has become cheaper, my impression is that nowadays, digital piano generation is mostly PCM samples recorded in multiple key (-zones) and multiple strike velocities, with a minimum of resynthesis to smooth things out. I might be wrong.
https://en.m.wikipedia.org/wiki/Inharmonicity


Are the frequencies a histogram in the frequency domain? I'm new to signal processing, and was learning about wavelets from this site.

In the graphs above, the time domain signal is obvious. The x axis shows the signal varying according to time. I don't understand the frequency domain though (even after watching a few videos).How does the frequency graph relate to the time graph above? Is it showing the frequency value for every amplitude value of a different signal (assuming "amplitude" is "Module" of the y axis) of the temporal graph?

Part of my assumption is obtained from this answer. But the question still remains, since if there are two signals of the same frequency but different amplitude, wouldn't they have a representation on the frequency graph where it's not possible to distinguish between them?

I'm even more confused because this example shows 20 millisecond samples being taken from a Fourier spectrum and being fed to a neural network like as though the x axis of the Fourier spectrum was the time domain rather than a "Module" value as shown in the graph above.

ps: I'd appreciate a plain-English explanation rather than one supplemented with equations that I won't understand.

Are the frequencies a histogram in the frequency domain?

No.

A histogram is the visualisation of a sequence of discrete counts. For example, a population pyramid can tell you the distribution of ages in a population.

The frequency spectrum represents strength. It doesn't have to be integer and it does not really mean that it represents counts of sinusoids at a specific frequency. The most trivial counter example of this is what is known as the Direct Current (DC) value which is frequency 0 (Zero) Hertz.

Although, you could possibly formulate a view of the frequency spectrum as some form of histogram (counts of elementary parcels of sinusoids), it would be very difficult to fit this view to the average value across the sample (which is what 0 Hz represents).

So, the frequency spectrum is not a histogram, it has nothing to do with a histogram. But what is it?

It is difficult to tell you without going through equations and properties because they represent our understanding of how the world works (for example, linearity and how things can be decomposed to parts and recomposed from parts). But here is an attempt:

This is a string trio:

Notice that the instruments have more or less the same shape but their size changes. Now, without going into equations, take it as a rule of thumb: The bigger the instrument the "lower" (in frequency) is its sound. (And, they can get...big).

When they all play together, they make up one continuous stream of sound that hits one's ears or a microphone.

The accurate analogy here is that the frequency spectrum of the sound recording would tell you how loud was each one of these 'violins' playing.

Joseph Fourier worked out a finer analogy that is very useful in DSP and goes by the name Discrete Fourier Transform (DFT).

He postulated that given any periodic waveform, there is a way to decompose it down to an infinite sum of sinusoids of varying strength and phase.

This is what the DFT does and from the DFT we derive various frequency spectra (namely, amplitude, power, phase).

  ...since if there are two signals of the same frequency but different amplitude, wouldn't they have a representation on the frequency graph where it's not possible to distinguish between them?

You see how this histogram view of the frequency spectrum doesn't really work now I hope. Two signals of the same frequency contribute to the strength of the same frequency in "the mix" and therefore do not appear as two individual measurements.

  I'm even more confused because this example shows 20 millisecond samples being taken from a Fourier spectrum...

But what you show is not a frequency spectrum. What you show is a spectrogram.

Without going into equations: The DFT decomposes a periodic waveform (any waveform) into an infinite sum of sinusoids. Because of the properties of those sinusoids, the DFT can tell you that a sinusoid is present in the mix but it cannot tell you WHEN it was present.

So, if you had a 1 hour recording of just two violins, a big one and a small one but with the small one just performing for 2 seconds, the DFT will still show high frequencies but it will not be able to inform you WHICH 2 SECOND SEGMENT of the 1 hour recording was it that the small violin was louder.

For this reason, the Short Time Fourier Transform (STFT) was developed. What this does is decompose the signal down to small (possibly) overlapping windows and obtain the DFT's of those windows.

This leads to an intermediate representation of the signal where we now can tell (within limits) which frequencies were present at which times.

This is what the last diagram in your question is depicting and this is why it still retains an axis of "Time". In that diagram, Time runs left to right and frequency runs bottom to top.

Hope this helps.

PS: I think that I understand why you don't feel like going through the equations but if you make a start, first with understanding the mathematical notation and then towards reading and understanding what it is that the equations are trying to describe, then I think that you will be greatly benefited in the long term.


get poles and zeros of frequency response I am working on a python based LTSPICE project.

I would like to get poles and zeros of AC simulation data. Is there a way to get them under use of the magnitude and phase out of the frequency response of my system?

This is a very complicated problems and I don't think there exists a one-size-fits-all solution. You can try Matlab's $invfreqz()$ and see if it works for your purposes https://www.mathworks.com/help/signal/ref/invfreqz.html

In general this is a error minimization problem but the actual data and the way you set the your error function and the search strategy makes a huge difference in the outcome. Things to watch out for

How good is your data? What's the frequency coverage range, density and shape of frequency grid, signal to noise ratio, etc.
What do you care about in the result. Dips vs peaks ? Phase vs amplitude, ?
Stability and causality constraints. How do you deal with non-causalities in the data (real or noise induced)
Search space representation, search strategy & algorithms.
Error surface shaping and optimization.
Convergence tracking, stop criteria, metrics for success (or lack thereof)

I developed one of these that is optimized around the specific constraints and requirements of my application and found invfreqz() to be totally useless. Your mileage may vary.


Numerical error while implementing a periodic pulse signal I'm trying to illustrate the principle of constructing a periodic signal using an elementary pattern. This is the code I use (python 2.7 + numpy + scipy.signal) :

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack as sfft
import scipy.signal as spsig

plt.close('all')

npts = 2**14 #number of points for both time-domain signals and FFT's
tsim = 100e-9 #ending time (in s)
fs = npts/tsim #sampling frequency

t = np.linspace(0,tsim,npts) #time bins vector
t1 = 30e-12 #pulse rise time
t2 = 80e-12 #pulse fall time
delta_t = 200e-12 #pulse width (including t1 and t2)
fsig = 1e9 #pulse repetition frequency
A = 1. #pulse amplitude

# here I construct the elementary pulse signal, it works fine
sig_elem = np.zeros(t.size)
sig_elem[(t>=0)&(t<t1)] = (A/t1)*t[(t>=0)&(t<t1)]
sig_elem[(t>=t1)&(t<delta_t-t2)] = A
sig_elem[(t>=delta_t-t2)&(t<delta_t)] = (-A/t2)*t[(t>=delta_t-t2)&(t<delta_t)] + (A/t2)*delta_t

periode = 1/fsig #this is the pulse repetition period
comb = np.reshape(np.array([np.round(t*1e9, 2)%1==0]), (npts,)) #this is where I construct the comb, I use numerical tricks to make it work
sig = spsig.fftconvolve(sig_elem, comb, 'full') #and the convolution itself, which should give a periodization of my elementary signal

##### FIGURE ####
plt.figure(1)
plt.plot(t*1e9, sig_elem, 'r',lw = 0.8)
plt.xlabel('t (ns)')
plt.xlim(-10,100)
plt.ylim(-1,5)
plt.grid(True)

# single sided FFT of the elementary signal. I restrict the plot frequency at 50 GHz since the sampling frequency is insanely high
sig_elem_fft = np.abs(sfft.fftshift(sfft.fft(sig_elem, npts)))
plt.figure(2)
plt.plot(sfft.fftshift(sfft.fftfreq(npts, 1/fs))*1e-9, sig_elem_fft, 'b',lw = 0.8)
plt.xlabel('f (GHz)')
plt.xlim(0,50)
plt.grid(True)

#plot of the dirac "comb". It seems good
plt.figure(3)
plt.scatter(t*1e9, comb)
plt.grid(True)
plt.xlim(-10,110)
plt.ylim(-0.5,1.5)

# plot of the convolve result, this is where I don't understand what's happening
plt.figure(4)
plt.plot(sig, 'r',lw = 0.8)
plt.xlabel('t (ns)')
#plt.xlim(-10,100)
#plt.ylim(-1,5)
plt.grid(True)

What I'm doing is defining a trapezoidal pulse :

$$
\begin{cases}
s(t) =
\frac{At}{t_1} & \text{ if }0 \leq t < t_1\\
A & \text{ if } t_1\leq t < \Delta t-t_2\\
\frac{A}{t_2}\left(\Delta t - t\right) & \text{ if } \Delta t-t_2 \leq t < \Delta t
\end{cases}
$$

$t_1$, $t_2$ and $\Delta t$ being respectively rise time, fall time and total pulse width. A is an arbitrary amplitude.
I implement it on $2^{14}$ (16384) points, making the total duration equal to 100 ns. Which makes the sampling frequency around 164 GHz.

This is what I get, the rest of the signal being equal to zero.

After creating this elementary pulse, I create a "dirac comb", using the same number of points, and using a 3 GHz repetition frequency.
Then, both signals are convolved using an FFT method (meaning both signals FFT's are calculated, multiplied element-wise and the result IFFT is computed afterwards).

What I should get is a periodic signal consisting of $s(t)$ defined above being repeated each 1 ns.
However, I get a strange signal, where the amplitude is sometimes A, sometimes 2*A. The period is however consistent.

Your "numeric trick" is not working. After running above code, have a look at :

plt.figure(); plt.plot(comb); plt.xlim((324,336)); plt.ylim((-0.1, 1.1))

Your "Dirac" train has peaks with width 2, hence you have the double-height amplitude pulses.

You can fix this by e.g. doing

comb2 = np.zeros(npts)
step = 1e-9*fs  # distance between peaks in terms of samples (floating point)
peakPos = step*np.arange(int(npts / step))  # integer position of the peaks
comb2[np.round(peakPos).astype(int)] = 1  # set the peaks
plt.figure()
plt.plot(comb)
plt.plot(comb2); plt.ylim((-0.1, 1.1)); plt.xlim((324, 336))


How can I mimic this ADC for a given resolution, sampling rate and input with MATLAB or Python? I'm quite new to the area so I need help to progress at some point.

Basically I have a signed 16 bit ADC which means it maps -10V to +10V range to 16 bit. So the max voltage +10V is 2^15-1 = 32767; zero is 0; and -10V is -2^15 = -32768.

I want to create/plot a sinusoid in MATLAB or Python where I can mimic exactly the ADC which samples a 1V 100Hz sinusoid input for 1 seconds which is sampled with 512Hz sampling rate. At the end I want to obtain the same FFT plot where I would obtain after a real ADC sampling.

ADCs are devices which can be simplified to behave ideally but are actually not doing so due to their physical limitations. Hence if you would like to make a realistic simulation of ADCs, then you should seriously consider their physical charactheristics as well. Otherwise if you are ok with the ideal mathematical definition of their outputs then the following simple MATLAB/OCTAVE code implements an example uniform quantizer. Please be very careful with the concept of simulation when it comes to ADC and quantizers.

clc; clear all; close all;

% ANALOG SIGNAL:
T = 1;          % observaton interval in seconds
f0 = 5;         % input sine frequency in Hz.
Fs = 1000;      % sampling frequency in Hz.

t = 0:1/Fs:T;   % sampling instants time

A = 1.25;       % amplitude of the input sine wave.
ax = A*cos(2*pi*f0*t)-(1E-15);   % ideal analog signal simulated with 64 bit IEEE floating point format
ax(ax <= -A) = -A + 1E-14;       % slight conditioning on the input so that decision limits are maintained.

figure,plot(t,ax);
title('simulated analog signal in 64-bit binary64 format');

% ADC :
N = 3;          % number of bits of uniform, symmetric, midtreat quantizer.
M = 2^N-1;      % number of levels: odd
qx = (2*A/M) *  round( (M/(2*A)) * ax ) ;   % QUANTIZED... (simulated within 64 bit container format)

figure,plot(t,qx),
title('simualted quantized signal');

The input and output are as follows:


Motivation for using gabor transform instead of fourier transform I am trying to understand the motivation for the above question for signal processing and image processing? Can anyone answer this question, or provide me with resources for this?

The fourier transform gives you very fine resolution in the frequency domain, but during the transformation, you loose all the information about when (for time signals) or where (for images) these frequencies occur in your input signal.

The Gabor transform alleviates this problem by windowing the base functions of the fourier transform with a Gaussian impulse. The transform gives the similarity between the signal and the time- and frequency-shifted window, so you may get a sense when/where certain frequencies occur in your input signal.

However, the windowing comes at the cost of losing some of the frequency resolution. The Gabor  transform is optimum in the sense that it's windowing function yields the minimum product between time and frequency resolution.


Evaluating regular Fourier Transform from DFT Suppose we know the DFT of a discrete limited sequence, some $X[k],\  k = 0, 1,\dots ,N-1$. How can we calculate the Fourier Transform of the same signal for a random frequency $\Omega$?

EDIT:

How about this, we can get the starting signal with inverse DFT, like: $$x[n] = \frac{1}{N}\sum\limits_{k=0}^{N-1} X[k]e^{j\frac{2\pi}{N}kn}$$
and then get the Fourier Transform normally with:
$$\begin{align}
X(e^{j\Omega}) &= \sum\limits_{n = -\infty}^{+\infty} x[n]e^{-j\Omega n}\\ &= \frac{1}{N}\sum\limits_{n = -\infty}^{+\infty} e^{-j\Omega n}\sum\limits_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}kn} \end{align}$$

Is this correct?

EDIT2:
As @Matt L. correctly states, I want the Discrete-time Fourier transform rather than the normal Fourier Transform, of course. In my native-language textbook Fourier Transform means both FT for continuous signals and DTFT for discrete signals. I hoped that that was clear from the edit. :)

Your question makes perfect sense, yet your solution will not work. Note that the expression for $x[n]$

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

is periodic with period $N$. This means that the sum

$$
\frac{1}{N}\sum\limits_{n = -\infty}^{+\infty} e^{-j\Omega n}\sum\limits_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}kn}$$

does not converge. Let's assume that $x[n]$ is defined for $0\le n< N$ and is zero otherwise. Then we get for the discrete-time Fourier transform of $x[n]$

$$X(e^{j\Omega})=\sum_{n=0}^{N-1}x[n]e^{-jn\Omega}=
\frac{1}{N}\sum_{n=0}^{N-1}\sum\limits_{k=0}^{N-1}X[k]e^{j\frac{2\pi}{N}kn}e^{-jn\Omega}\tag{1}$$

This formula actually works and can be used if only the DFT $X[k]$ is given. Of course this only works for finite length signals $x[n]$. Equation (1) simply means that the DFT contains the complete information about the (finite-length) signal $x[n]$ and any frequency point $\Omega$ of the DTFT can be evaluated using the discrete values $X[k]$.


Separability vs. Resolution - Synonymous terms? In terms of discrete representation of a variable, let's say time measurements in absolute values in [seconds], is there a difference between separability and resolution? When do speak of either aspect? Or are those terms used synonymously?

  Example: Measuring arrival time of two physically distinct flights
  landing, are there examplary different contexts, in which we would
  distinguish between time separability of flight arrival times and time
  resolution?

Most technical terms, such as these two, do not get their definitions purely from their etymology, but rather from the context of application, by experience and by tradition of acceptance.

And for this case, your understanding of separability as the allowence (or ability) of someone to distinguish between two closest time of arrivals seems a synonym to the more classical term resolution in your context. Nevertheless, they will (and can) take their absolute precise definitions only when they are given mathematical expressions to be quantified...

Note that, in multi-dimensional signal processing, the term separability assumes a completely different meaning too, and which is not related to the discussion here...


Carrier recovery samples Should carrier recovery be done always 1 sample per symbol? This is ok if carrier recovery is before timing. But Nco inside carrier recovery may work at higher rate

If my carrier recovery is first then should i do at n samples/symbol. ie 4
If number of samples / symbol are more is there any advantage.

This depends on your modulation choice and approach to carrier recovery. But consider this for a specific case of a raised cosine pulse shaped signal:

Observe (using the eye diagram would be simple) the significant symbol to symbol jitter everywhere except in close proximity to the symbol decisions — for this case I would opt to use 1 sample per symbol (the decisions) to estimate symbol to symbol phase since those samples have the least jitter (this would be the samples after all equalization, timing recovery and matched filtering is complete, just prior to decision where the complex phase can still be measured between two successive symbols).

This approach works well if you are using the Gardner TED for symbol clock recovery since it can function well with relatively large carrier offsets.

You can evaluate your specific waveform with the same considerations to determine what approach would have the minimum estimation error (and not the error on any given sample but the rms error within your carrier tracking loop bandwidth).

Here is an example shown below for a raised-cosine filtered waveform. Typically the filtering is split between the transmitter and receiver with Root-Raised Cosine (RRC) filters such that sufficient pulse shaping for reducing out of band emissions is provided while allowing for a matched filter for optimum SNR (against white noise conditions) in the receiver. In the plots below you see the signal as received (if noise-free) and then after it is passed through the RRC matched filter.  Notice the wider variation in the zero crossings after the RRC filter in the plot on the right compared to the plot on the left. It is the zero crossing locations that are used, on average, to determine the data clock. This is why, without additional processing to minimize the zero crossing variation (which can be done), it is better to use the the waveform on the left for timing recover. However notice the sharp points at the correct sampling locations for symbol decision on the plot on the right compared to the one on the right-- these are the best samples (at one sample per symbol using the best derived sample from timing recovery) to measure phase rotation for carrier recovery given the Gardner TED can operate with significant carrier offset.


How do I implement a footstep recognition algorithm? One that recognizes the start and end points of each footstep? I have some time series data captured from a person's footsteps/strides (specifically, a person on rollerblades). It came from an IMU sensor placed on a person's boot. Each data point captures acceleration, angular velocity, etc.
Here's a graph of "acceleration in the x direction vs. time" for one foot:

Problem/Goal: I want to write an algorithm to detect all the footsteps in the data, specifically, the start and end times of each footstep. (Start = foot leaves the ground; End = foot touches the ground).
Here's an example below. Each shaded region represents the duration of a footstep (from start to end). Annotated data like this is available for me to use, perhaps in a supervised learning model.

Here are some of the approaches I am considering after searching around. (Feel free to SKIP this section.)

Dynamic sliding window + classifier

This involves somehow sliding a dynamic window over the time series. At each window, run a classification model to determine whether the windowed data is a footstep or not. A dynamic window is necessary because each step varies in length.
From this blog, I know how to create a signal classification model. I can use Fast Fourier Transform or Wavelet Transform to extract features to use in a supervised learning model.

Gradient Boosting Machine

I read this paper that created a model to recognize the start/end points of footsteps...
for turkeys lol. It uses a Gradient Boosting Machine and trained on data where all data points belonging to a "footstep" were annotated (e.g. all shaded areas in the graph above).

Question: Does anyone have any suggestions for how to approach this problem? Or perhaps any comments on whether the approaches above make sense or don't make sense.
Note: My knowledge in Deep Learning is pretty basic, so I am not considering Deep Learning approaches for the sake of time. Perhaps something for the future.

This signal looks pretty clean.  I would stick to time domain feature extraction.  Maybe 1st and 2nd derivative of time domain data.  If you're just counting, you can almost just use the magnitude and see that these steps have the most negative values or acceleration.  You could use the negative peaks to center a search for the edges of the "step", and to discriminate against the other negative peaks that are in between the steps.


A basic question about the use of moving average vs low-pass filters in DSP I'm newbie novice in digital signal processing. If one has a signal with some noise superimposed in time series, for which type of information/analysis one would use moving average or LP filter in practice.

Imagine you have a 5 seconds of sampled data in time where the sampling rate is 200Hz.

Can you give a practical example when would you use LP filter and when would you use a moving average filter? I read that moving average is a type of LP filter but I would like to see the differences in big picture.

I see many engineers use MATLAB or Python functions and perform Low Pass filter using FFT but not moving average. So this made me to ask this question.

A moving average filter can be thought of as a type of low-pass filter that doesn't have any control over its bandwidth for a fixed number of taps.

For a finite impulse response (FIR) filter, the output signal $y[n]$ is given in terms of the input signal $x[n]$ and the filter taps $h[n]$:
$$
y[k] = \sum_{n=0}^{N-1}h[n]x[k-n].
$$
The filter length in this case is $N$ "taps" (or $N$ coefficients). A moving average filter has coefficients that are all equal:
$$
h[n] = \frac{1}{N}, \qquad n = 0, 1, \ldots, N-1,
$$
whereas in general, a low-pass filter (LPF), can have different values for each tap. This allows you to control the frequency selectivity of the filter.

Here is an example showing the difference, for a system with a sample rate of 200 Hz and a filter with 32 taps:

The LPF was designed to have a cutoff frequency of 50 Hz, but it can be set to any value you want. With a moving average filter the filter is narrowly focused around the 0 Hz component ("DC"), and the peak gets narrower the more taps you have in the filter. Another problem with using a moving average filter as an LPF is that it has high sidelobes (the ripples to either side of the main peak) compared to a "properly designed" filter.

The taps for the LPF are:

[-0.00116417 -0.00139094  0.00195951  0.00293134 -0.00437535 -0.00637313
  0.00902803  0.01248116 -0.0169409  -0.02273977  0.03045372  0.04118039
 -0.05729369 -0.08500841  0.14720004  0.45005217  0.45005217  0.14720004
 -0.08500841 -0.05729369  0.04118039  0.03045372 -0.02273977 -0.0169409
  0.01248116  0.00902803 -0.00637313 -0.00437535  0.00293134  0.00195951
 -0.00139094 -0.00116417]

Python code for generating this figure:

from pylab import *
from scipy.signal import windows, firwin

Fs = 200.

N = 32
h_lpf = firwin(N, 50, nyq=Fs/2., window='hamming')
h_ma = ones(N)*1./N

print h_lpf

M = 512

X_lpf = fftshift(abs(fft(h_lpf, M)))
X_lpf /= X_lpf.max()
X_ma = fftshift(abs(fft(h_ma, M)))
X_ma /= X_ma.max()

f = arange(M)*Fs/M - Fs/2.

figure()
plot(f, X_lpf, f, X_ma)
xlabel('$f$ (Hz)')
ylabel('$|H(f)|$')
legend(('LPF', 'Moving Average'))
grid()
show()


Figuring magnitude and phase response I've got a linear time-invariant system $$y[n]=\frac{8}{9}y[n-1]+x[n]$$
which I transformed into a transfer function $$Y(z)=\frac{8}{9}Y(z)*z^{-1}+X(z) =>\frac{Y(z)}{X(z)}=\frac{1}{1-\frac{8}{9}*z^{-1}}=\frac{z}{z-\frac{8}{9}}$$
The pole is positive, which would make the system stabile.
However, I am having some trouble getting the magnitude/phase responses $|H(\omega)|$ and $\theta(\omega)$.
How should I proceed with the system to figure these out? I've read about using the substitution for $z$ as$$Y(z)=\frac{z}{z-\frac{8}{9}}=\frac{e^{iw}}{e^{iw}-\frac{8}{9}}$$
or even continue to transform the exponentials to their $\cos$ + $\sin$ forms, but i am not sure how i should proceed with them.

As an addition, how would the system change if I would to swap $x[n]$ with a sin or a cosine of some value? I assume it would be the same path towards the responses but with different values.

$$
|H(z)| = \sqrt{\Re^2 H(z) + \Im^2 H(z)} \\
\theta \big(H(z)\big) = \arctan \frac{\Im H(z)}{\Re H(z)} \\
$$

$\Re z$ is the real part of z;

$\Im z$ is imaginary part.

For a single-pole filter $y_n = b_0 x_n - a_1 y_{n-1}$,

$$
H(z) = \frac{b_0}{1 + a_1 z^{-1}}; \\
\Re H(z) = b_0 \frac{1 + a_1 \cos \omega}{1 + 2 a_1 \cos \omega + a_1^2}; \\
\Im H(z) = b_0 \frac{a_1 \sin \omega}{1 + 2 a_1 \cos \omega + a_1^2}.
$$

With $b_0 = 1$ and $a_1 = -\frac{8}{9}$, your filter may be unstable.

For a low-pass filter, try $b_0 = \frac{1}{9}, a_1 = -\frac{8}{9}$.


Least Squares Linear Phase FIR Filter Design In explication ''the geometric interpretation of least squares'' Typically, the number of frequency constraints is much greater than the number of design variables (filter coefficients). In these cases, we have an overdetermined system of equations (more equations than unknowns). Therefore, we cannot generally satisfy all the equations, and are left with minimizing some error criterion to find the ``optimal compromise'' solution.'

What is meant by ''constraints''
What is the relation between constraints and the number of design variables (filter coefficients).

In this case must be using a least-squares approximation.But, are there other cases when constraints less than or equal to the number of design variables?

I think that the term "constraints" is not a very fortunate choice in this context, but what is meant is the number of frequency points that are specified:

$$H(\omega_i)\stackrel{!}{=}D(\omega_i),\qquad i=1,2,\ldots,K\tag{1}$$

where $H(\omega_i)$ is the actual frequency response evaluated at frequency $\omega_i$, $D(\omega_i)$ is the desired response at $\omega_i$, and $\stackrel{!}{=}$ means that this is not an actual equality but we're trying to minimize some norm of the error $H(\omega_i)-D(\omega_i)$ over all frequencies $\omega_i$, $i=1,2,\ldots,K$.

Generally, the number of "equations" in $(1)$ is much greater than the number of filter coefficients $N$. In practice one often choses $K$ $5$ to $10$ times greater than $N$. This is why you get an overdetermined system of equations which can be solved in a least-squares (or any other) sense.

If $K<N$ you have an interpolation problem, which means that normally all equations can be satisfied with equality, but you have no control over what happens in between the specified frequency points $\omega_i$. Also, if $K<N$ the filter coefficients are not completely specified because you have an underdetermined system of equations. Such a system can be solved with an additional requirement such as that the vector of filter coefficients has minimum norm. But traditionally, filter design problems are not formulated as underdetermined systems of linear equations.


Educational purpose - What is the correct way to simulate a multipath fading channel which has ISI I am following the codes for rayleigh fading channel given in the webpage: http://www.raymaps.com/index.php/m-qam-bit-error-rate-in-rayleigh-fading/

I am confused and need help in clarifying some basic questions regarding how to simulate a Rayleigh fading channel with $L$ taps? This is not clearly mentioned.

Question 1:  I have studied that a flat fading channel is one which has 1 tap. In the webpage, a flat fading is simulated by h=(1/sqrt(2))*(randn(N,1)+1i*randn(N,1)); where $N$ denotes the number of bits or the number of data points. So, there are $N$  taps and definitely this is not a 1 tap channel. Please correct me where wrong.

Question 2: Also, how to simulate a fading and flat fading with say $L = 5$ taps?

This is how  I have done but I am not sure how to properly do such that there is inter-symbol-interference as well. If number of paths, $L$ through which the signal reaches the receiver with some delay among the different paths. If the delay spread is greater than the symbol duration then we get inter symbol interference. In the following code, I am not sure if inter-symbol-interference is there or not and if I am correclty applying the theory to simulate the fading channel.

clear all
% Number of data points
N = 100;
L= 5; % number of taps
%Data symbols of 0/1
x = randn(N,1)>0;
%channel coefficients as Rayleigh random variable with variance 0.5 for real and imaginary component
h = 0.5*(randn(L,1)+1i*randn(L,1));
%Channel is an FIR filter
y = filter(h,1,x);

Lets say we want to transmit a sequence of discrete data $\left\lbrace x[n] \right\rbrace$. But because we are living in analog world, the sequence must be modulated.

Call $T_s$ is symbol duration and use a set of orthonormal waveforms $\left\lbrace p_n(t) = p(t-nT_s), n \in \mathbb{Z} \right\rbrace$, (baseband) signal $x(t)$ can be written as
\begin{equation}
x(t) = \sum_{n} x[n] p(t - n T_s) \tag 1
\end{equation}

The equation (1) can be thought of modulating $x[n]$ to the n-th dimension of $x(t)$. To get back $x[n]$, we project $x(t)$ to this n-th dimension.

This projection is implemented by passing the received signal over the matched filter $q(t) = p^*(-t)$ and then sampling at $t=n T_s$. This projection is perfect if $g(t) = p(t) \star q(t)$ satisfies the the criterion $g(t) = 1$ if $t=0$ and $g(t) = 0$ if $t=m T_s\textrm{, } m \in \mathbb{Z}\textrm{, } m \neq 0$, i.e. Nyquist ISI criterion. Otherwise, the project will introduce error and by projecting to the n-th dimension, we don't get only $x[n]$ but also $x[m \neq n]$.

($\star$ means convolution)

This is bandlimited-caused ISI, not what you are asking but the model is necessary to understand the second ISI type, multipath-propagation-caused ISI.

Another source of ISI is multipath phenomenon. The received signal can be written as

\begin{equation}
y(t) = \sum_{i} a_i(t)x(t-\tau_i(t)) + w(t) \tag 2
\end{equation}
where $i$ is the index of physical path, $a_i(t)$ is its complex gain and $\tau_i(t)$ is the delay introduced by the path. $w(t)$ is AWG noise at receiver.

Apply the matched filter and sampling the output at $t = m T_s$,
\begin{eqnarray}
y[m] = y(t) \star q(t)|_{t = m T_s} = \sum_{i} a_i(m T_s) \sum_{n} x[n] g(m T_s - n T_s - \tau_i(m T_s)) + w[m] \\
y[m] \stackrel{l=m-n}{=} \sum_{l}x[m-l]h_l[m] + w[m]
\end{eqnarray}

with the channel tap $h_l[m]$, i.e. channel tap $l$ at the sample time $t=mT_s$,

\begin{equation}
h_l[m] = \sum_{i}a_i(m T_s)\times g(l T_s - \tau_i(m T_s)) \tag 3
\end{equation}

Now assume that the channel impulse reponse is invariant (a relaxed assumption could be that channel taps does not change during a coherence time, we call it underspread assumption),

$$y[m] = \sum_{l}x[m-l] \times h_l + w[m] \tag 4$$

The equation (4) is what you call FIR model. The equation (3) is the model of channel tap.

Remember Nyquist ISI criterion above ? $g(t) = 1$ if $t=0$ and $g(t) = 0$ if $t=m T_s\textrm{, } m \in \mathbb{Z}\textrm{, } m \neq 0$. Thus the power of the composite filter $g(t)$ must "concentrate" between $[-T_s/2, T_s/2]$ as in image for raised cosine pulse.

Thus the term $g(l T_s - \tau_i(m T_s))$ in Equation (3) is significant if and only if $-T_s/2 < |l \times T_s - \tau_i(m T_s)| < +T_s/2$.

  We can interpret as the channel tap $l$ is contributed by the physical
  path $i$ whose delay $\tau_i$ is around $t = l \times T_s$.

  The channel impulse reponse can be thought being sampled by symbol duration $T_s$. And the number of tap depends on the relation between symbol duration $T_s$ and delay spread $\tau_m = \mathrm{max}_i (\tau_i)$, i.e. $L \approx \tau_m / T_s$.

  As soon as $L > 1$, $y[m]$ is contributed by several $x[m], x[m-1], ...$ and it is multipath ISI. An equalizer will be needed to mitigate the effect of this ISI.

  If $\forall \tau_i \ll T_s$, all physical path gains $a_i$ will contribute to only the channel tap $l=0$, we have flat fading, or single tap model.

What you are doing with matlab is playing with Equation (4) and you don't need to care about symbol duration and delay spread anymore.


Fourier transform before calculating the phase When and why do we need to perform Fourier transform before calculating the phase of the image?

The discrete Fourier transform (DFT) does by definition transform your image into the frequency domain, mapping it to a matrix of complex numbers.

Again, by definition, a complex number has a phase and a magnitude component.

Thus, you don't "need to perform FT before calculating the phase of an image", you get the phase of the image as a result of the FT.

It's a bit like saying "why do I need to perform an addition before calculating the sum of 2 number?".

I would suggest you to learn more about Fourier Transform


Cross Correlation in Java (Android) Using FFT (Frequency Domain Convolution / Correlation) I want to implement cross-correlation with the start-symbol and the signal using FFT.

I zero-pad both signals to length N = x.len + y.len -1
Using the convolution theorem:

        corr(x, y) = IFFT( FFT(x_padded) .* conj(FFT(y_padded)))

Do you have any good library recommendations I could use to compute both the FFT and IFFT?
I am running this on an Android app and efficiency is important.

Moreover, the best correlation is the highest peak in the array of the result. How does the index k of the array corr(x,y) relate to x? Does that mean the best alignment is when y lags k samples behind x?

If index k refer to the highest correlation coeffitient and let n be the length of each x and y, so would n-k be the lag for y at maxima of the correlation function.

For some applications I prefer the cosine similarity instead of plain correlation, but which is less efficient in calculation.


Filter Coefficients in MATLAB are always float data type I will be using the filter coefficients generated by MATLAB in my code that will run on a microcontroller. The poor microcontroller is very bad in floating point Arithmetic. So What do I do? Shall I directly truncate all float values to nearest int? shall I do ceil? or floor? Is there any function that can automatically do this conversion for me?

Update
All my filter coefficeints are coming like: 0.465, 0.76, 0.23 etc, so If I use ciel then all will become zero. What should I do now?

Actual Filter Coefficients are:

> -0.00385638
>  0.004944457
>  0.01505063
>  0.018768283
>  0.009635631
> -0.01192891
> -0.035809426
> -0.045043857
> -0.023993426
>  0.032106934
>  0.111777547
>  0.18997355
>  0.238374966
>  0.238374966
>  0.18997355
>  0.111777547
>  0.032106934
> -0.023993426
> -0.045043857
> -0.035809426
> -0.01192891
> 0.009635631
> 0.018768283
> 0.01505063
> 0.004944457
> -0.00385638

Normally you would use fixed-point arithmetic for performing DSP on a a floating-point-challenged CPU.

So for, say, 16 bit fixed-point coefficients with a range of -1.0 to +1.0 your coeffs would translate to e.g.:

0.465    =>    0.465 * 0x7fff = 0x3b5

If your microcontroller does not have direct support for fixed point arithemtic (most DSPs do, general purpose microcontrollers typically do not) then you'll need to take care of scaling when you multiply or divide. (Addition and subtraction work as normal of course, provided you are not mixing different fixed point types.)


solution on the time domain becomes “periodic” after the inverse fourier transform I was trying to solve european option pricing problem using Conv method (introduced by Lord in 2008 https://pdfs.semanticscholar.org/0632/460bd50b2151f74ac40028df4cc60e73a884.pdf).

The final step of this algorithm is to transform the solution on the frequency domain to the time domain by inverse fast Fourier transform, which means that I can retrieve approximations for all the grid points.

However, I noticed that only the solutions in the center of the time domain (in this case, x=40, and I used log grid points) are accurate. Here is the plot: the orange line represents the approximation using FFT, the blue line represents the real solution.

Moreover, the approximation looks periodic. So I was wondering what kind of situation will inverse FFT give a 'periodic-like' solution? I did some research on this topic, and foundthat DFT assumes the input is periodic, which may cause  "time domain aliasing" (https://www.dspguide.com/ch10/3.htm).

I also found that "spectral leakage" may be related to my problem. Can anyone give me some advice on this? Thank you very much.

P.S. By periodic, I mean the solutions on the two ends are very close to each other. I changed lots of parameters, and this always holds.

Multiplying a kernel and signal spectrum in Fourier domain lead to a circular convolution and not a linear convolution, so in order to it become linear convolution you must zero pad your signal and kernel before taking the Fourier transform (up to M+N-1 where M is the signal's length and N is the kernels's length). (if you compare the blue and orange signals, the beginning of orange signal has higher values because of of higher value at the end of blue signal and it has lower values in the end because of lower values in the beginning of the blue signal).


Expected value of signals If I have two signals $x$ and $y$ and this expression $E\{x[n]y[n-k]\}$, what is the expected value of the product?

I have a basic understanding of EV based on probability, such as $$E[X]=\sum_{i=1}^{\inf}x_ip_i$$ and I understand that EV of flipping coin is 0.5 if $x_1=0$ and $x_2=1$ and both probabilities are 0.5 or for rolling a dice is 3.5 etc.

However, when I see the expression like this the one above, I find it confusing. I don't know for instance where the probability term is.

This is the cross-correlation function of the two discrete-time random processes $x[n]$ and $y[n]$. It measures the similarity of the processes for a given time lag between the two. If $x[n]$ and $y[n]$ are jointly wide-sense stationary, then the cross-correlation only depends on the lag $k$, and not on the absolute time index $n$.

In terms of probabilities (or rather probability density functions), the expectation $E\{x[n]y[n-k]\}$ is given by

$$E\{x[n]y[n-k]\}=\int_{-\infty}^{\infty}xyf_{XY}(x,y;n,n-k)dxdy$$

where $f_{XY}(x,y;n_1,n_2)$ is the joint probability density function of $x[n_1]$ and $y[n_2]$.

In practice, we usually hope that the processes are ergodic, which allows us to estimate such an expectation by taking the appropriate time averages.


How to quantify delay in filtered output? I have a signal containing two sine tones 30 kHz and 60 KHz. I used a second order IIR notch filter to filter out the 60 KHz component from the signal, so now the signal has only 30 KHz component.

When comparing the filtered out signal with a pure 30 KHz tone generated, the earlier samples are very much different but later samples are matching to some good extent.

Plot with earlier samples:

Plot with later samples:

What are the reasons why there is a delay before response being proper?
How can i mathematically quantify this?

Like after how many samples the filtered output will be reasonably matching with that of pure tone?

It's the transient and the group delay which are associated with the initial conditions and the LTI filter phase response $\phi(w)$ respectively.

The group delay (in samples) associated with the LTI filter is:
$$\tau = - \frac{d \phi(w) }{dw} $$

The output will be shifted by $\tau$ samples wrt to the input. You shall compansate this shift for synchronization.

This delay is constant for linear phase filters which are of symmetric FIR type filters. For IIR type filters, since their phase response is not linear therefore getting a fixed group delay is not possible.


Reconstruction of a signal from non-uniform samples So I'm having some troubles with a signal that's been sampled but where the samples don't have consistent intervals. I've been looking into the NFFT (in python) as a starting point but I'm a complete noob with signal processing so I'm not really getting anything out of it.

So I have one array with the sample values and another with the corresponding timestamps. But I really have no idea how to get the spectrum out of it. I've done the first part as described here.

Sorry about the sparse question.

The problem you have to solve is having a model for how your samples uniquely represent actually continuous input.  For equidistant samples, the usual model is to stipulate (and then guarantee by analog filtering, oversampling, digital filtering, and decimation) that there are no signal components exceeding the Nyquist frequency.

What is your model going to be for the kind of non-uniform sampling you are working with?  If the information you work with is not sufficient for that model, you are going to have non-uniquely interpretable results, namely aliasing.


Spectral efficiency of with OFDM vs without OFDM? OFDM has advantage over frequency selective channel with addition of cyclic prefix after IFFT of MPSK or MQAM symbol. but I wonder that the spectral efficiency of with OFDM vs without OFDM.

The spectral efficiency of OFDM is strictly worse than that of a pulse-shaped QAM signal with the same rate.

OFDM requires a guard interval, on which no useful information is transmitted.

It is common to dedicate several subcarriers to pilot signals.

It is often not pulse-shaped (or, more accurately, pulse-shaped with a rectangular pulse).

Since each non-pilot subcarrier is essentially a narrowband QAM signal, you can see that its spectral efficiency cannot be better. However, the point of OFDM is to allow for low-complexity frequency-domain equalization, with spectral efficiency given a less important role.