Fourier filtering in Matlab

Hi to everyone, I am using MATLAB for my research in Atmospheric Science. I want to filter my daily data (wind data) using band pass filter of 20-90 days. My problem is that i have only 90 days which is not enough to pass the 20-90 days variabilities. I used Foufilter.m by putting samplingtime=1, center frequency = 11/360 (i calculated center frequency as (1/20+1/90)/2,frequency width=7/360 (1/20-1/90),shape=20,and mode=0 (for band pass). Also i pad the data with 365 zeros on each side. so i have 820 days in total including actual data. But i am getting straight line. I do not know where problem lies? Can anyone suggest anything? Also i want to control the wave number in the filter, e.g. i want to retain the wave number 1-5 in the data (I am working on tropical area and wave number 1-5 is MJO oscillation).

[reddy]

 Hello all and Waqar;

I am also using MATLAB for Ocean Science research. As Waqar told he want to filter 20-90 day wind data, similarly I am having daily data of current for 117 days(everyday one measurement). Now I want to see POWER spectrum after 20 days means (mask or filter below 20 day signals). For that I just used Wayne King's code. But I the plot still shows peaks with in 20 days where I dont want to see them with in 20 days. I want them to be flattened. Hope you all got my point. And I am attaching the sample data set as well as the code, output image.

Please check my code and let me know where did I went wrong.

u = xlsread('sample.xlsx');
d = fdesign.lowpass('N,Fp,Ap,Ast',20,1/20,1,20,1);
D = design(d);

 y = filter(D,u);
N= length(u);        % number of points
T = N;
t=(0:N-1)/N; % defining time
f = (0:N/2-1)/T;        
per = 1./f;         % finding the corresponding frequency and period in days
p_sur_cur= abs(fft(y))/(N/2);           % absolute value of the FFT
pow_sur = p_sur_cur(1:(N)/2).^2;        % take the power of positive freq half
plot(per,log(pow_sur),'color','m','LineWidth',2); legend('U-current','Location','southeast');grid on

As we can see the peaks in the 0-20 days, which can not be shown if we do filter. That is where I lost...

Thanks in advance

Answer 1

[Wayne King]

Unfortunately padding the time series with zeros is not going to improve your ability to recover oscillations for which your original time series is insufficiently long. Why don't you just lowpass filter your data for with passband period of 20 days? As you correctly note, one problem is that even for an oscillation of 20 days, you only enough data for slightly more than 4 cycles. Do you have the Signal Processing Toolbox? If so, you can use fdesign.lowpass. I'll assume from your post that your sampling interval is 1 sample/day.

 D = fdesign.lowpass('N,Fp,Ap,Ast',20,1/20,1,20,1);
 d = design(D);

Then to filter your data:

output = filter(d,input);

I'll just simulate some data for you:

 t = 0:89;
 x = cos(2*pi*1/40*t)+cos(2*pi*1/20*t)+randn(size(t));
 y = filter(d,x);
 plot(t,x,'b'); hold on;
 plot(t,y,'r','linewidth',2);
 % for zero-phase filtering
 y1 = filtfilt(d.Numerator,1,x);
 plot(t,y1,'k-.','linewidth',2)

To look at the spectrum:

 ydft = fft(y);
 ydft = ydft(1:length(y)/2+1);
 freq = 0:1/90:1/2;
 figure;
 plot(freq,abs(ydft))
 xlabel('Cycles/day'); ylabel('Magnitude')

[Waqar]

Hi Thanks for your kind response. I need to use band-pass filter for my time series. I cannot use low-pass filter. Also, i consulted this issue with my friend and we came to the understanding that i need to reconstruct the time series using wave number 1-5. Here wave number 1 corresponds to 90 days whereas the wave number 5 corresponds to 18 days. Do you know how to reconstruct the time series in MATLAB using wave number 1-5? One possibility is to use fft but i do not know how to use my time series while applying filter?

Answer 2

[Waqar]

Hi Thanks for your kind response. I need to use band-pass filter for my time series. I cannot use low-pass filter. Also, i consulted this issue with my friend and we came to the understanding that i need to reconstruct the time series using wave number 1-5. Here wave number 1 corresponds to 90 days whereas the wave number 5 corresponds to 18 days. Do you know how to reconstruct the time series in MATLAB using wave number 1-5? One possibility is to use fft but i do not know how to use my time series while applying filter?