Selecting a part of the frequency spectrum

4 次查看（过去 30 天）
giacomo labbri
giacomo labbri 2020-7-28
评论： Star Strider 2020-7-28
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');

请先登录，再回答此问题。

回答（1 个）

Star Strider
Star Strider 2020-7-28
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 个评论
giacomo labbri
giacomo labbri 2020-7-28
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 2020-7-28
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by