power spectral density of nonstationary data

I am hoping those with seasoned experience with signaling can assist me with two questions. The ultimate goal is to obtain a power spectral density of very nonstationay optical data sampled at 1 kHz (attached). There is a huge DC offset that persists even after detrending or subtracting the mean. For such nonstationary data, pwelch seems most appropriate but a rather odd looking series of deminishing regularly space peaks is outputted across the frequency domain. Is the windowing wrong? Is this the wrong analysis to do?

First, here is the pwelch based code I ran:

T = readtable('time_series_data.csv');
time_series = T.time_series;
x = time_series - mean(time_series);
fs = 1000;
[pxx_smooth, f] = pwelch(x, hann(length(x)/2), [], fs);
figure;
plot(f, 10*log10(pxx_smooth), 'LineWidth', 1.5);

If you run the code on the attached time-series you get no distinct peaks and a huge off wtih a very narrow frequency domain. This makes no sense.

Second, despite supposedly not correct to do since this is non-stationary data, I ran an FFT:

fs = 1000;
%x = detrend(time_series);
x = time_series-mean(time_series);
figure; 
L = length(x); 
plot(fs/L*(-L/2:L/2-1),abs(fftshift(x)),"LineWidth",2)

With an FFT you get more defined peaks and over a much larger frequency range (see attached jpg), but confused why the asymmetry when using ffshift, where it should be symmetrical about x = 0. Why the asymmetry?

I am convinced I doing something fundamentally wrong with the processing or windowing and odd looking output is not inherently due to the data itself. Any assistance in how to correctly attain a PSD for this data and what is the most valid approach for nonstationary data of this nature would be greatly appreciated.

Answer 1

[Paul]

Hi hxen,

When uploading the data in a .csv file, please also include the code that reads the data from that file. It wasn't hard to figure out in this case, but it makes it easier for anyone who wants to help.

The second code snippet didn't actually take the FFT of the data. See below for correction.

Do the plots now look as expected?

T = readtable('time_series_data.csv');
time_series = T.time_series;

x = time_series - mean(time_series);
fs = 1000;
[pxx_smooth, f] = pwelch(x, hann(length(x)/2), [], fs);
figure;
plot(f, 10*log10(pxx_smooth), 'LineWidth', 1.5);

ans = 9.3922e+03

fs = 1000;
%x = detrend(time_series);
x = time_series-mean(time_series);
figure; 
L = length(x);
%plot(fs/L*(-L/2:L/2-1),abs(fftshift(x)),"LineWidth",2)
% should the data be converted to dB?
plot(fs/L*(-L/2:L/2-1),abs(fftshift(fft(x))),"LineWidth",2)

copyobj(gca,figure);
xlim([-3,3])

[hxen]

Hello Paul. Thanks for pointing out the inadvertant code ommission for the reading in the data. I appreciate you replying, but this doesn't address what I had posed: the non-stationarity.

The codes you have are identical to what I had submitted so unfortunately not addressing the issue. Nonstationary data is expected to yield a dominant peak about zero for an FFT and as mentioned despite detrending or subtracting the mean, this persists for both an FFT and welch, which is particularly designed to deal wtih non-stationary data and addess the edge effect. Yet clearly that is not happening.

In addition, I had to truncate the data to the <5MB allotted upload limit. The shorter recording time-series I provided actually exaccerbates the non-stationarity in the ouputs. But there are should be frequency ranges that are dominant outside of zero. My suspicion as first mentioned is the windowing for (Hann or Hamming) is not chosen correctly for this data type and applied correctly.

Edit: Also, I should add the limtied frequecy output for the welch analsys is also an issue.There should be frequency ranges with peaks out well over 20 Hz.

[Paul]

I might come back to this later, but want to point out that the code I provided is NOT identical to the code in the question. The code in the question did not actually call fft in the second method, which is why the plot looked (and was) incorrect.

Can you explain or provide a reference to what you mean by "non-stationarity." I just want to make sure we are on the same page.

[hxen]

Hi Paul. The data is non-ergodic. Any statistical measures is essentially random over time, including ampltitude. There are frequent long segments of no activity that is no doubt accounting for the dominant peak about zero. Welch is suppose to handle these kinds of cases and why I am convinced doing the windowing wrong. Han or Hamming windowing supposedly dampens these edge effects. It ulitmately be too the data is just inherently not ideal for this kind of analysis. Thanks for your time.

[Paul]

Full disclosure: I'm a bit out of my comfort zone when it comes to random processes and spectral estimation. Make of that what you will.

Sounds like "statistical measures ... random over time" describes the underlying random process from which this data is sampled as being non-stationary. "The data is non-ergodic." I thought ergodic means something else, though, IIRC, one might be a sufficient condition to imply the other. It's been a while.

Anyway, lets look at the data:

T = readtable('time_series_data.csv');
time_series = T.time_series;
x = time_series - mean(time_series);
fs = 1000;
plot((0:numel(x)-1)/fs,x),grid

It certainly doesn't look the same to the eye across different windows, so I guess that suggests the underlying random process is not stationary.

"Welch is suppose to handle these kinds of [non-stationary] cases" To the contrary, I thought Welch's method makes the fundmental assumption that the data is sampled from a wide-sense, stationary random process as stated at pwelch: "Because the process is wide-sense stationary ..." Nevertheless, I think Welch's method can be used if the process isn't exactly stationary, though closeness to stationary might be open to interpretation (see disclosure above).

Here's the code using pwelch, but with a correction.

Referring to the linked doc page, in the four-argument form as used in the original code

[pxx_smooth, f] = pwelch(x, hann(length(x)/2), [], fs);

the fourth argument is nfft. So fs (1000) is being used as nfft, and the output frequency vector, f, is normalized frequency (identified as w on the doc page) spanning 0-pi (as shown above)

The correct call (assuming defaults for noverlap and nfft) is:

[pxx_smooth, f] = pwelch(x, hann(length(x)/2), [], [], fs);
figure;
plot(f, 10*log10(pxx_smooth), 'LineWidth', 1.5);

●

Now the frequency axis goes to fs/2, as it should.

Here's the one-sided fft aproach, with the fft squared. I did not scale to convert to a psd for direct comparison to pxx_smooth.

figure; 
L = length(x);
X = fft(x);
X2 = X.*conj(X);
plot(fs/L*(0:L/2),10*log10(X2(1:L/2+1)),"LineWidth",2);ylim([-300,150]);

Of course dB(X2) is not defined at f = 0 because mean(x) == 0 (effectively). However, pxx_smooth will be finite at f = 0 because the windowed segments do not all have zero mean.

The output of pwelch can vary drastically depending on the input parameters. Compare the output when using the default window, which as I understand is a Hamming window of length that obtains eight segments with the default 50% overlap, whereas the pxx_smooth used the Hann window with three segments with 50% overlap (I think). I don't know why the psd estimates diverge from each other for f > 10 Hz (see disclosure).

[pxx2, f2] = pwelch(x, [], [], [], fs);
figure;
plot(f, 10*log10(pxx_smooth), f2, 10*log10(pxx2), 'LineWidth', 1.5);grid
xline(10);

[hxen]

Thank you and appreciate your time Paul.

[Paul]

You're quite welcome. I'm find this all very interesting to explore.

Read in the data

T = readtable('time_series_data.csv');
time_series = T.time_series;
x = time_series - mean(time_series);
fs = 1000;

Call pwelch with three different windows and compare the results.

[P{1},f] = pwelch(x, hann(length(x)/2), [], [], fs);
 P{2}    = pwelch(x, hamming(length(x)/2), [], [], fs);
 P{3}    = pwelch(x, rectwin(length(x)/2), [], [], fs);
figure;
plot(f,db(horzcat(P{:}))),grid,legend('hann','hamming','rectwin');

copyobj(gca,figure);
xlim([190,210]);

For reasons I don't understand, the hann window yields a much different result than the others. However, the hamming and rectwin also both yield the sinc-like characteristic.

Divide the signal into three segments of length N/2, with 50% overlap. I think these would be the segments on which the pwelch calls above operates and averages the results.

N = numel(x);
xx{1} = x(1:N/2);
xx{2} = x(N/4+1:3*N/4);
xx{3} = x(N/2+1:end);

Plot transform of each window and the power spectrum of each windowed segment.

figure
t = tiledlayout(3,5);
win = {@hann,@hamming,@rectwin};
cellfun(@(win) helper(xx,win,fs),win,'Uni',false);

Consistent with pwelch, the power spectrum with the hann window is not like the others, even though its spectrum appears, by eye, to be very similar to that of hamming. If anything, one might expect the rectwin to be the outlier. Uncler why the results with hann are so different than with hamming.

Let's try pwelch, but with much shorter windows.

[P{1},f] = pwelch(x, hann(4096), [], [], fs);
 P{2}    = pwelch(x, hamming(4096), [], [], fs);
 P{3}    = pwelch(x, rectwin(4096), [], [], fs);
figure;
plot(f,db(horzcat(P{:}))),grid,legend('hann','hamming','rectwin');

I guess this plot shows the effect of a lot more averaging, but still unclear to me why hann results in so much more attenuation for f > 10 Hz.

The fact that hann and hamming look very similar (to me) in the frequency domain, but have much different effects on the results, is very interesting. I mean, multiplication in the time domain is convolution in the frequency domain, but convolving the same segment with very similar window transforms yields much different results. I don't know if that's an effect in general, or unique to this particular input signal.

function helper(xx,win,fs)
w = win(numel(xx{1}));
W = fft(w);
N = numel(W);
f = (-N/2:N/2-1)/N*fs;
xx = cellfun(@(xx) xx.*w,xx,'Uni',false);
nexttile;plot(f,db(abs(fftshift(W)))),grid;ylabel('abs'),title(functions(win).function);
nexttile;plot(f,angle(fftshift(W))),grid,ylabel('angle')
XX = cellfun(@(xx) fft(xx),xx,'Uni',false);
XX = cellfun(@(XX) XX.*conj(XX),XX,'Uni',false);
cellfun(@(XX) plot(nexttile,f,db(fftshift(XX))),XX,'Uni',false);
end

[hxen]

Just seeing this. Been reading about the selection of the window spacing—a critical consideration. Rec is the worse for this kind of data and Hann more ideal as it attentuates the edges. However the apparent attenuation > 10Hz range is due to the spectral leakage <10 Hz due to many closely spaced freqs wtih high amplitudes there. Your visualization of the different windowing is very helfpul. Hamming is better for freq resolution but Hann better for amplitude extraction. Agree too this is a fascinating subject worth investing more time into learning. I've been fortunate until now that regular FFT analysis has worked perfectly fine. Again, appreciate you and your time!

[Paul]

The frequency domain plots of the window functions I showed above are misleading. Zero padding is needed to reveal their actual structure.

N = 100;
whann    = hann(N);
whamming = hamming(N);
f = (-N/2:N/2-1).'/N;
Whann    = fft(whann);
Whamming = fft(whamming);

N5 = N*5;
f5 = (-N5/2:N5/2-1).'/N5;
Whann5    = fft(whann,N5);
Whamming5 = fft(whamming,N5);

figure
plot(f,db(abs(fftshift([Whann,Whamming],1))),f5,db(abs(fftshift([Whann5,Whamming5],1)))),grid
ylim([-125,50]);
legend('Whann','Whamming','Whann-padded','Whamming-padded')