Why do I see a drop in the last datapoint (Nyquist frequency) of the spectra derived from pwelch?

1 view (last 30 days)
Just to provide some background, I have been working with velocity data and using the pwelch function to obtain the spectra I need. Generally, I would work with datasets with sampling frequencies of 4 or 8 Hz. This time I am working with a slightly different dataset, sampled at 16 Hz with a different technology. I have been getting a drop in the last datapoint of the spectra, which relates to the Nyquist frequency. I see that in individual spectra but it becomes even more obvious when I average together spectra derived when the average current velocity was very slow, as in the example:
Here is a piece of the code I have been using, which actually uses the pwelch function:
nfft = 256;
window = 'hann';
win = hann(nfft);
overlap = 0.5;
noverlap = overlap * nfft;
tmp_vel_X = ADV_timetable.vel_X(i1:i2);
tmp_vel_X = detrend(tmp_vel_X);
[ADV_Processed.spectra_X(:,i_burst), ADV_Processed.f] = pwelch(tmp_vel_X, win, noverlap, nfft, CFG.fs);
It is in a loop so that the whole signal is cropped between i1 and i2. Each cropped section has 9600 datapoints and CFG.fs = 16.
I'm not very experienced with signal processing methods, I have used only a few things. So I believe maybe I'm missing out on something that could be very basic for someone who masters the topic. I have been looking up other forums and signal processing books but I couldn't find something that helps with that specific question.
Unfortunately I'm not allowed to share a sample of the data...
Thanks a lot in advance to anyone who tries to give me some light! :)
  2 Comments
Mathieu NOE
Mathieu NOE on 1 Mar 2023
hello
as we cannot use your data , I simply put together this simple demo using the same parameters as yours
I created a sufficiently long signal so that we get a significant amount of averaging by pwelch
the signal itself is a combination of a 1 Hz tone plus some broadband noise (random)
you see in this case the baseline spctrum is flat up to the Nyquist frequency (Fs/2 = 8 Hz here) , so my believe is that if you see a drop in your signals periodogram , it's because your signal simply does not have any energy at this frequency.
maybe this was hiden before when you were sampling at lower rate (8 or 4 Hz).
Fs = 16;
t = 0:1/Fs:500;
x = sin(2*pi*t)+randn(size(t));
nfft = 256;
win = hann(nfft);
overlap = 0.5;
noverlap = overlap * nfft;
pwelch(x, win, noverlap, nfft, Fs)
Larissa Perez
Larissa Perez on 1 Mar 2023
Hi Mathieu,
Thanks a lot for taking the time to look into this. I appreciate it.
I'm not sure if that's it. The point of this analysis is looking at turbulence, which actually represents itself in many frequencies, higher than 8Hz.
Yesterday I went through some old work of mine and I found the same sudden drop in the spectra, even with lower sampling frequencies. I just had never noticed it before because it wasn't so clear as it is in this case.
I wonder if it is some sort of artifact or if maybe I'm just not using the optimal parameters for pwelch in my case. I have tried using different nfft but it just makes it less visible in the plot. When i look at the values, the drop at 8Hz is still there.

Sign in to comment.

Accepted Answer

MarKf
MarKf on 1 Mar 2023
It may be because it's the Nyquist frequency, rolloff/aliasing/edge artefacts should happen afterwards at higher freqs but in some cases at Nyquist depending on the system. Also because the A/D converter of the new data acquisition system could have an anti-aliasing low pass filter, sometimes they have a cutoff frequency just before (like 90%) the Nyquist frequency

More Answers (1)

Bruno Luong
Bruno Luong on 1 Mar 2023
Edited: Bruno Luong on 2 Mar 2023
No I beg to differ the answers that have given to you. The reason is the convention of onesided spectrum. In the code computepsd.m (called by pwelch) the Nyquist and DC are halft of other frequency by convention of one-sided FFT
The original code (for even nfft as with your case) is
Sxx = [Sxx_unscaled(1,:); 2*Sxx_unscaled(2:end-1,:); Sxx_unscaled(end,:)]; %
If you replace this statement with (what I called "non-standard one-sided FFT")
Sxx = [2*Sxx_unscaled(1,:); 2*Sxx_unscaled(2:end-1,:); 2*Sxx_unscaled(end,:)];
I beg you don't see any drop with this mofified code.
As example I show the example based on @Mathieu NOE
Fs = 16;
t = 0:1/Fs:5001;
x = sin(pi*t);
x(end) = [];
nfft = 256;
win = hann(nfft);
overlap = 0.5;
noverlap = overlap * nfft;
pwelch(x, win, noverlap, nfft, Fs)
%Sxx = [2*Sxx_unscaled(1,:); 2*Sxx_unscaled(2:end-1,:); 2*Sxx_unscaled(end,:)]
The scaling by factor of 2 of the no-extreme part of the spectrum raises the question of correct "interpretation" that is always subject of many errors if users blindly apply blackbox code and don't pay enough attention (or ignore) to such detail.
  5 Comments
Bruno Luong
Bruno Luong on 2 Mar 2023
Edited: Bruno Luong on 2 Mar 2023
Note also that if I trick pwelch by convert x to complex; then the two-sided spectrum is returned without the central part boosted. So the amplitude of Nyquest frequency looks continuous to the rest.
Fs = 16;
t = 0:1/Fs:5001;
x = sin(pi*t);
x(end) = [];
nfft = 256;
win = hann(nfft);
overlap = 0.5;
noverlap = round(overlap * nfft);
pwelch(complex(x), win, noverlap, nfft, Fs); % <= trick MATLAB here
xlim([0 Fs/2])
MarKf
MarKf on 2 Mar 2023
That makes sense and it is a thing that OP should definitey try with their data. I dismissingly considered this as "edge artefacts", which is a lazy way to ignore stuff that happens with empirical data, but yes, they are usually caused of the way PSD is calculated on discrete/finite signals.
But the sharp cut off seen in the figure is not just a roll off but looks like a bandpass, by my experience this is usually caused by the system itself.

Sign in to comment.

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!