在MATLAB中重建FFT的時間序列

我想重建來自FFT的太陽黑子信號，時間序列和週期圖在以下站點http://www.mathworks.com/help/matlab/examples/using-fft.html。我寫了下面的代碼，但結果並沒有類似於原波：在MATLAB中重建FFT的時間序列

YY=Y(1:floor(n/2)) 
% magnitude 
mag_fft = 2*abs(YY)/length(Y); 
% phase angle 
ang_fft = angle(YY); 


[new_mag,new_i]=sort(mag_fft,'descend'); 
new_ang=ang_fft(new_i); 
new_freq=freq(new_i) 
wave=zeros(1,length(YY)); 
wave=new_mag(1); 
t=1:length(YY) 

for(i=1:70) 
wave=wave+new_mag(i).*sin(2*pi*new_freq(i)*t+new_ang(i)); 
end 
wave=wave-mag_fft(1) 
figure;plot(year(t),wave,'-b') 
hold on;plot(year(t),relNums(t),'-r')

什麼想法？

來源

2014-02-23 kernel

在命令行輸入'help ifft'。 –

謝謝，是的，ifft做fft的倒數。我想我想要做的是一樣的，但可能是不同的方法。我試圖計算每種模式的幅度和相位，然後再次組合它們以獲得原始時間序列。 – kernel

%http://www.mathworks.com/help/matlab/examples/using-fft.html 
% sunspots 
% sunspots have period of 10 years 
%% 
clc;clear all;close all; 
load sunspot.dat 
year=sunspot(:,1); 
relNums=sunspot(:,2); 
figure;plot(year,relNums) 
title('Sunspot Data') 
plot(year(1:50),relNums(1:50),'b.-'); 

yfft = fft(relNums);%figure;plot(ifft(yfft)-data1d,'r') 
%yfft = fft(data1d); iyfft=ifft(yfft); 
[sum(relNums) yfft(1)] 
yfft(1)=[]; % we grid rid of the first value as it corresponeding to zero frequency. 
N=length(yfft)+1; 
yfft=yfft.*2./N; 
%% 
power_fft = abs(yfft);power1_fft = sqrt(yfft.*conj(yfft)); 
figure;plot(power_fft,'-b');hold on;plot(power_fft,'rO') 

ang_fft = angle(yfft);real_fft= real(yfft);imag_fft= imag(yfft); 
figure;plot(real_fft);hold on;plot(imag_fft,'-r') 
figure;plot(angle(yfft)) 
    ph = (180/pi)*unwrap(ang_fft); % phase in degrees 

% Now the total length of the per and all other powers should be N-1 because there is no 
% more corresponding poweres and phases, and the number of frequencies before the nequiest is 
    Nneq=length(N./(1:N/2)); 

    Nm1=N-1; per=N./(1:Nm1); freq=1./per; 
    [per'/12 power_fft(1:Nm1)/100 ] % so as to display the period in years 

    %% ytyt 
    ndat=length(relNums); 
    x=0:ndat-1; 
    sumharmony1(1:Nneq,1:ndat)=0; 
    sumharmony2(1:Nneq,1:ndat)=0; 

    for i=1:Nneq 
    % those two forms are equal, the last one is called the cos form. 
    %  sumharmony1(i,:)=sumharmony1(i,:)+real_fft(i)*cos(2*pi*x/(per(i)))- imag_fft(i)*sin(2*pi*x/(per(i))); 
    sumharmony1(i,:)=sumharmony1(i,:)+power_fft(i)*cos(2*pi*x./(per(i))+ang_fft(i)); 

    end 
    y1 =sum(relNums)/N+ sum(sumharmony1); 
    %y2 =sum(tmp)/N+ sum(sumharmony2); 

    figure;plot(relNums);hold on; plot(y1, 'r'); 
figure;plot((relNums-y1')) % However, the excellent results, we couldnot yet reach to the that of the built in function ifft. 
    figure;plot(relNums(1:100),'-ob');hold on; plot(y1(1:100), 'r'); 
    % note that we multiply by 2 because of using the window hanning.enter code here

來源

2016-12-12 06:23:26 kernel

在MATLAB中重建FFT的時間序列

回答

相關問題