Issues with power spectral density calculation of linear ramp

Hi, I am trying to use some old code in a new application and running into a strange issue. My code takes a 1D array of data, splits it into N samples, and averages the PSD (actually the linear spectral desnity - sqrt of PSD) of each to get a smoothed value.

My custom code uses only fft(), it's basically an implementation of https://mathworks.com/help/signal/ug/power-spectral-density-estimates-using-fft.html

When I used to run this code (~2014), the averaging worked as expected. When I ran it today on new data, there was an offset on low-frequency data that's associated with a linear ramp. This is a problem because I use the amplitude at low frequency a part of interpretation. I haven't changed the code.

On investigating, I've simplified as much as possible and can show the same effect with the matlab function 'periodogram'. To illustrate the problem we don't need to stack/smooth- we just use the first 1/10th of the data, and we can already see the offset that's bothering me. Full code below.

My understanding is that the power spectral density of a portion of a continuous signal should be the same as the complete signal- there should be no offset anywere, let alone one that impacts part of the signal (low-frequency ramp) compared to another part (the white noise portion at high frequency).

Any advice is appreciated.

% random noise + linear ramp with N samples
N = 1000000;
y = randn(1,N)'+(10/N)*(1:N)';
% take 1/10th of y as a sample
y_10 = y(1:N/10);
% sampling is at 1s so don't need to create time array

figure(1);
clf;
hold on;
ylabel('Pressure (psi)')
xlabel('Time (seconds)')
plot(y);
plot(y_10);
legend 'data' '1/10th of data' 

figure(2);
clf;
hold on;
grid on;
xlabel('Frequency (Hz)')
ylabel('Noise spectral density (psi/rt.Hz)')
set(gca, 'YScale', 'log')
set(gca, 'XScale', 'log')

% I'm plotting linear spectral density, but you see the same thing with psd
[pxx,f] = periodogram(y);
plot(f,sqrt(pxx));
[pxx,f] = periodogram(y_10);
plot(f,sqrt(pxx));
legend 'From data' 'From 1/10th of data'

Plots:

[Jason MILNE]

To give a concrete example of why this might be a problem- in electronic noise there is concept of 'corner frequency' where 1/f noise (ramp) takes over from gaussian noise (flat). If I tried to extract the corner frequency from this data, I would get different results for the two curves. I'm using a linear ramp rather than 1/f noise so it's not exactly the same, but it illustrates the issue.

Answer 1

[dpb]

y=detrend(y);