MATLAB Answers

Selecting a part of the frequency spectrum

1 view (last 30 days)
giacomo labbri
giacomo labbri on 28 Jul 2020
Commented: Star Strider on 28 Jul 2020
Hi,
I am trying to characterize a daily cycle of some variable (temperature etc...). To do that I am taking the timeseries converting it into frequency domain (using fft) and then selecting only the frequncies greater or equal than the daily frequency. Then I am converting back to space domain the selected frequencies (using ifft) and I would expect to get the daily cycle and all the shorter variation (for example the variations over an hour). But I don't see the daily cycle, I just see the shorter variations. To make my problem simpler I tested my code on a test function (sum on three sines: one with a period of 5 days one with a period of 1 day and one with a period of half on hour). I am attaching my code
My main question is: am I doing something wrong in the code?
Thanks in advance,
Giacomo
t=1:46080; %minutes
y = sin(2*pi *t /30 )+sin(2*pi *t / (24*60) )+sin(2*pi *t / (5*24*60)); %test function
y=transpose(y);
Fs =1/60; % Sampling Frequency
Fn = Fs/2; % Nyquist Frequency
L = numel(t); % Signal Lengths
ym=y-mean(y); % Subtract Column Means From All Columns (Eliminates D-C Offset Effect)
FTy = fft(ym); % Fourier Transform (if you want it Scaled For Length uncomment the division by L)
Fv = linspace(0, 1, fix(L/2))*Fn; % Frequency Vector
Iv = 1:numel(Fv); % Index Vector
%Y=ifft(FTy);
Fv_dc=Fv(Fv>=(1/(24*3600))); %selectin the frequencies greater or equal then the daily freq.
Iv_dc=Iv(Fv>=(1/(24*3600)));
FTy_dc=FTy(Iv_dc);
dc=ifft(FTy_dc,L); %this should be the dirunal cycle and all the shorter period variation
dc=real(dc);
sp=(abs(FTy(Iv))).*2; %full power spectrum
sp_dc=(abs(FTy_dc)).*2; %power spectrum of the frequencies greater or equal to the diurnal frequency
figure(1)
plot(t,y,t,dc) %original series and the series correponding to the frequqncies greater or equal to the daily freq.
figure(2)
loglog(Fv,sp,'.-b',Fv_dc,sp_dc,'.-r') % full spectrum and selected spectrum
xline(1/(24*3600),'--r');

  0 Comments

Sign in to comment.

Answers (1)

Star Strider
Star Strider on 28 Jul 2020
It is difficult to follow what you are doing, however it appears that you are ignoring half of the fft results, and that is causing the problems that you see. (The fft produces symmetrical results, with one half being the complex conjugate of the other half. You can see this most easily by using fftshift and then plotting the result as the function of an appropriate, symmetrical frequency vector.)
Forget about filtering in the frequency domain. Just use the highpass or bandpass functions (depending on what you want to do) to filter in the time domain, and you will likely get the result you want. (These functions are in the Signal Processing Toolbox, R2018a and later versions.) If you have an earlier version of the Toolbox, designing filters in MATLAB is straightforward nusing the available funcitons.
.

  2 Comments

giacomo labbri
giacomo labbri on 28 Jul 2020
Thanks for your imput. I tried to use the highpass function but the result I get are not strange. I asked a separate question on the highpass function but I have recived no answer yet ( https://it.mathworks.com/matlabcentral/answers/562322-highpass-function-vs-manual-filtering?s_tid=prof_contriblnk )
If you have any advice on how to correct/make more redable my code they would be very much welcome : )
Star Strider
Star Strider on 28 Jul 2020
My pleasure!
I have no idea what you want to do, so you need to design the fillters to produce the result you want. Use the fft result to guide your choice of the highpass passband frequency. Filter design usually requires a bit of experimentation.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!