%傅立葉法
%本方法的作用是：對指定信號進行頻譜分析
%
%   輸入參數：singal 為待分析的信號；輸入應為256個值的數據 
%            當ISI給-1時，做時間抽選FFT正變換；當ISI給1時，做頻率抽選FFT正變換
%   輸出參數：H 為目標函數幅值；f 為頻率，
%


%%%%%%%%%%%%%%%  傅立葉變換  %%%%%%%%%%%%%%%%%
filein='D:\data_signal';
ISI = -1;
FID = fopen('D:\data_signal.txt','r')
    [ N ] = fscanf(FID,'%d\n',1);
    [f,count] = fscanf(FID,'%f  %f\n',[2,N]);
    f=f';
fclose(FID);  

% %%%%%%%%%%%%%%  反傅立葉變換  %%%%%%%%%%%%%%%%%
% filein='D:\data_spectrum';
% ISI = 1;
% FID = fopen('D:\data_spectrum.txt','r')
%     [ N ] = fscanf(FID,'%d\n',1);
%     [f,count] = fscanf(FID,'%f  %f\n',[2,N]);
%     f=f';
% fclose(FID); 

%%%%%%%%%%%%%%%  基本參數設置  %%%%%%%%%%%%%%%%%
MMAX=1;
NN=1;
J=1;

%%%%%%%%%%%%%%%%  判斷數據長度是否是2的整冪數  %%%%%%%%%%%%%%%%
for I=1:15
    M=I;
    NN=NN*2;
    if (NN==N)
        break;
    end
end

%%%%%%%%%%%%%%%%  得到輸入數據的實部和虛部  %%%%%%%%%%%%%%%%
  for I=1:N
      data1(I)=f(I,1);                     %數據的實部
      data2(I)=f(I,2);                     %數據的虛部
      if ISI == -1
          signal(I)=data1(I);  %輸入數據的模
      end
      if ISI == 1
         signal(I)=sqrt(data1(I).^2+data2(I).^2);  %輸入數據的模
      end
  end
  
%%%%%%%%%%%%%%%%  數據的重排序  %%%%%%%%%%%%%%%%
for I=1:N
    if I<=J
        temp1=data1(J);
        temp2=data2(J);
        data1(J)=data1(I);
        data2(J)=data2(I);
        data1(I)=temp1;
        data2(I)=temp2;
    end
        M=N/2;    
        while  J>M  
            J=J-M;
            M=fix((M+1)/2);           
        end
            J=J+M;       
end

%%%%%%%%%%%%%%%%  主體計算  %%%%%%%%%%%%%%%%
while MMAX<N
    ISTEP=2*MMAX;
    for M=1:MMAX
        theta=pi*(ISI*(M-1))/MMAX;
        W1=cos(theta);
        W2=sin(theta);
        for I=M:ISTEP:N
            J=I+MMAX;
            temp1=W1*data1(J)-W2*data2(J);
            temp2=W1*data2(J)+W2*data1(J);
            data1(J)=data1(I)-temp1;
            data2(J)=data2(I)-temp2;
            data1(I)=data1(I)+temp1;
            data2(I)=data2(I)+temp2;
        end
    end
    MMAX=ISTEP;
end

%%%%%%%%%%%%%%%%  反傅立葉變換時  %%%%%%%%%%%%%%%%
if ISI>0
    for I=1:N
        data1(I)=data1(I)/N;
        data2(I)=data2(I)/N;
    end    
end

%%%%%%%%%  將正變換后的數據輸出為文本格式  %%%%%%%%%%%%%%
if ISI == -1
 FID = fopen('D:\data_spectrum.txt','w')
    fprintf(FID,'%20d \n',N);
    for I=1:N
        fprintf(FID,'%20.10f  %20.10f\n',data1(I),data2(I));
    end 
 fclose(FID);  
end

 %%%%%%%%%  將反變換后的數據輸出為文本格式  %%%%%%%%%%%%%%
if ISI == 1
  FID = fopen('D:\data_signal.txt','w')
     fprintf(FID,'%20d \n',N);
     for I=1:N
         fprintf(FID,'%20.10f  %20.10f\n',data1(I),data2(I));
     end 
  fclose(FID); 
end 

%%%%%%%%%%%%%%%%  求出振幅和相位  %%%%%%%%%%%%%%%%
for I=1:N
    H(I)=complex(data1(I),data2(I));           %  振幅的計算 
    module(I)=sqrt(data1(I)^2+data2(I)^2);     %  模的計算
    P(I)=atan(data2(I)/data1(I));              %  相位的計算
end
%%%%%%%%%%%%%%%%%%%%%%  畫圖  %%%%%%%%%%%%%%%%%%%%%%
if ISI == -1
    figure(1);
    subplot(311); %圖形顯示分割窗口 
    plot(signal);
    title('輸入信號');
    subplot(312); %圖形顯示分割窗口 
    plot(module);
    title('頻譜分析振幅譜');
    subplot(313); %圖形顯示分割窗口 
    plot(P);
    title('頻譜分析相位譜');
end
if ISI == 1
   figure(2);
   subplot(211); %圖形顯示分割窗口 
   plot(signal);
   title('輸入信號');
   subplot(212); %圖形顯示分割窗口 
   plot(module);
   title('反傅立葉變換后振幅譜');     
end