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

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! :)

[Mathieu NOE]

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]

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.

Answer 1

[MarKf]

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

[Larissa Perez]

I'll investigate that but it makes sense. Many thanks!!

Answer 2

[Bruno Luong]

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.

[Bruno Luong]

Note that the drop of last onesided frequancy (EDIT: not Nyquist as Paul's pointed out) is not occured for odd value of nfft

Fs = 16;
t = 0:1/Fs:5001;
x = sin(pi*t);
x(end) = [];
nfft        = 257;
win         = hann(nfft);
overlap     = 0.5;
noverlap    = round(overlap * nfft);
pwelch(x, win, noverlap, nfft, Fs)

[Paul]

Following up on this, the reason it doesn't drop off at Nyquist for odd value of nfft is that the plot actually stops before the Nyquist frequency when nfft is odd (and I asusme the equation that does the scaling of Sxx must be different as well) becasue the Nyquist frequency is not a bin when nfft is odd

Fs = 16;
t = 0:1/Fs:5001;
x = sin(pi*t);
x(end) = [];
nfft        = 257;
win         = hann(nfft);
overlap     = 0.5;
noverlap    = round(overlap * nfft);
pwelch(x, win, noverlap, nfft, Fs);
xlim([7.5 8])

[Bruno Luong]

@Paul, correct I should write last onesided frequency instead of Nyquist. I fix my comment, thanks.

[Bruno Luong]

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]

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.