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

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.

Answer 1

[Brian Neiswander]

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:

http://www.mathworks.com/help/dsp/ref/butter.html

http://www.mathworks.com/help/signal/ref/filtfilt.html

Hope this helps.

-Brian