Power Spectral density of function with an offset

I am trying to compute the power spectral density of 2 functions. They are similar in terms of the magnitude of their fluctuations from their mean, but one has a large offset. This leads me to believe that the function with the large offset should have a power spectral density above the function without an offset. However this is not what I see and I am hoping someone can explain this to me. I have attached my code which approximates the signals and then calculates the power spectral densities using several different methods. Methods 1 and 2 show that the power spectral densities or periodograms of the two functinos are on similar orders of magnitude, contrary to what I expect. Method 3 then attempts to use the autocorrelation function to compute the power spectrum but this doesnt work well at all. Finally if I use matlabs in built power spectrum function I see something like what I expect, that the function with an offset has a much higher power spectrum.

Can anyone explain to me how to reconcile these four methods of computing the power spectrum? I believe they should all return the same. I also would like to understand why the first two methods produce similar power spectrums for the two functions, when one clearly has much more power, owing to the offset.

Finally is there is a physical meaning if I take the fft of the square of a time series?

clear all 
close all
n=100000;
tmax=1e5;
t=linspace(0,tmax,n);
dt=mean(diff(t));
fs=1/dt;

S1=1e-5+5e-8*sin(2*pi*100*t)+randn(1,n)*1e-8;
S2=randn(1,n)*1e-8;
% 
% plot(S1)
% 
% figure(2)
% plot(S2)
% return
 
freqvec=linspace(0,fs,n/2);
 
%fft
ftS1=fft(S1);
ftS1=ftS1.*conj(ftS1).*(1/(fs*n));
ftS2=fft(S2);
ftS2=ftS2.*conj(ftS2).*(1/(fs*n));
 
figure(1)
loglog(freqvec,abs(ftS1(1:end/2)))
hold on
loglog(freqvec,abs(ftS2(1:end/2)))
 
 
%periodgram
[p,fv]=periodogram(S1,[],length(S1),fs);
[p2,fv2]=periodogram(S2,[],length(S2),fs);
figure(2)
loglog(fv,sqrt(p))
hold on
loglog(fv2,sqrt(p2))
 
 
%autocorrelation
AC1=autocorr(S1,'NumLags',length(S1)-2);
AC2=autocorr(S2,'NumLags',length(S2)-2);
ft1=fft(AC1);
ft2=fft(AC2);
figure(3)
freq=linspace(0,fs,n/2);
loglog(freq,abs(ft1(1:n/2)))
hold on
loglog(freq,abs(ft2(1:n/2)))
 
%power spectral density
freqsamp=mean(diff(freq));
spec1=pspectrum(S1,freqsamp);
spec2=pspectrum(S2,freqsamp);
 
figure(4)
loglog(spec1)
hold on
loglog(spec2)

Answer 1

[Mathieu NOE]

Hi

I modified your code. there were several bugs there.... but FYI an offset is just a DC value, so in the FFT (or psd) spectrum it will only afect the dc frequency bin , but not the higher bins (which are related to the dynamic content of your signal)

This is validated by the code. the difference in the plots are only for the DC frequency.

I commented the "autocorr" method section as I have not exactly the same version of matlab.

I think the 3 others computations are enough for the demo

All the best

n=1e6;      % samples
fs = 1000; % sampling freq (Hz)
signal-freq = 100; % make sure fs > 2.56 * signal freq

t=(0:n-1)/fs;
dt=1/fs;

tmax=max(t);

S2=randn(1,n)*1e-8;
S1=1e-5+5e-8*sin(2*pi*100*t)+S2;    % so S1 is based on S2 + other

% 
% plot(S1)
% 
% figure(2)
% plot(S2)
% return
 
freqvec=linspace(0,fs/2,n/2);
 
%fft
ftS1=fft(S1);
ftS1=ftS1.*conj(ftS1).*(1/(fs*n));
ftS2=fft(S2);
ftS2=ftS2.*conj(ftS2).*(1/(fs*n));
 
figure(1)
% loglog(freqvec,abs(ftS1(1:end/2)))
% hold on
% loglog(freqvec,abs(ftS2(1:end/2)))
semilogy(freqvec,abs(ftS1(1:end/2)),'b',freqvec,abs(ftS2(1:end/2)),'r')
legend('S1','S2');
 
%periodgram
[p,fv]=periodogram(S1,[],length(S1),fs);
[p2,fv2]=periodogram(S2,[],length(S2),fs);
figure(2)
% loglog(fv,sqrt(p))
% hold on
% loglog(fv2,sqrt(p2))
semilogy(fv,sqrt(p),'b',fv2,sqrt(p2),'r')
 legend('S1','S2');

% %autocorrelation
% AC1=autocorr(S1,'NumLags',length(S1)-2);
% AC2=autocorr(S2,'NumLags',length(S2)-2);
% ft1=fft(AC1);
% ft2=fft(AC2);
% figure(3)
% freq=linspace(0,fs,n/2);
% loglog(freq,abs(ft1(1:n/2)))
% hold on
% loglog(freq,abs(ft2(1:n/2)))

%  
%power spectral density
% freq=linspace(0,fs,n/2);
% freqsamp=mean(diff(freq));
% spec1=pspectrum(S1,freqsamp);
% spec2=pspectrum(S2,freqsamp);

nfft = 640;
%   [Pxx,F] = PWELCH(X,WINDOW,NOVERLAP,NFFT,Fs) returns a PSD computed as
%   a function of physical frequency (Hz).  Fs is the sampling frequency
%   specified in Hz.

[spec1,f]=pwelch(S1,hanning(nfft),0.75*nfft,nfft,fs);
[spec2,f]=pwelch(S2,hanning(nfft),0.75*nfft,nfft,fs);
figure(4)
semilogy(f,spec1,'b',f,spec2,'r')
legend('S1','S2');

[Lorcan Conlon]

okay, thank you for your response. I understand what you mean that the difference will be at DC. I have attached a figure below of the two temporal signals. Looking at them briefly can you guess whether you would expect them to have similar PSDs?

[Mathieu NOE]

first , I don't get the same plot as you

and , no , you cannot guess what the dfference on the psd can be from the time plot (IMHO)

when I am in doubt, I do time / frequency analysis to compare two signals. this way you can easily see in which time span / frequency range there could be differences.