re-build time series from pwelch after removing peaks in the spectra

5 views (last 30 days)
Richard
Richard on 19 Aug 2015
Answered: Brian Neiswander on 21 Aug 2015
I have time series data that is similar to the following:
t = 3650; % number of days in ten years
t = 1:t; % generate time vector
fs = 1; % sampling frequency (days)
A = 2; % amplitude
P = 365; % period (days), the time it takes for the signal to repeat itself
f1 = 1/P; % number of cycles per second (i.e. how often the signal repeats itself every second).
y = A*sin(2*pi*f1*t); % signal
y = awgn(y,15,'measured');
figure(1);
subplot(211);
plot(t,y)
% spectra for data
n = 1024;
Nfft = n;
[Pxx, Freq] = pwelch(y,hamming(n),n/2,Nfft,1);
subplot(212);
plot(log10(Freq),Freq.*abs(Pxx));
where I have 10 years of data that is characterised by a strong seasonal cycle, which is obvious from the data but also identified using pwelch. Next, I would like to re-generate the original series from the spectra but set the peaks equal to the power of a smooth background model (I'm assuming this would be the same as removing the peaks from the spectra and smoothing with the nearest points?), which would remove the seasonal cycle from the data and return the remaining time series. Note that the reason I'm doing this is because I will eventually have a more complicated series that will have additional periods and I would like to remove them aswell.

Answers (1)

Brian Neiswander
Brian Neiswander on 21 Aug 2015
Hi Richard,
I recommend using a band-stop filter rather than manually removing a peak from the spectra and then recreating the time-domain signal. Manually removing peaks can produce strange results when working with measurement data due to content bandwidths and harmonics.
In your example data above, you can use the "butter" function to design a Butterworth band-stop filter with cutoff frequencies at 0.8 to 1.2 times the seasonal cycle frequency:
>> [b,a] = butter(2,[0.8 1.2]*f1/(fs/2),'stop'); %2nd order band-stop filter
Then filter the data using the "filtfilt" function:
yf = filtfilt(b,a,y);
You can then check that the seasonal cycle has been attenuated in the power spectra:
>> [Pxxf, Freqf] = pwelch(yf,hamming(n),n/2,Nfft,1); %estimate psd of filtered signal
>> figure;
>> semilogx(Freq,Freq.*abs(Pxx));
>> hold on;
>> semilogx(Freqf,Freqf.*abs(Pxxf),'r');
If you have the Signal Processing Toolbox, you can also try using the interactive "Filter Design & Analysis Tool" to quickly test out different filters:
>> fdatool
You can follow the links below for more information on the "butter" and "filtfilt" functions:
Hope this helps.
-Brian

Community Treasure Hunt

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

Start Hunting!