MATLAB-Fourier變換

 

•  將振幅為1Hz的正弦波和振幅為0.5的5Hz正弦波相加后進行Fourier分析，研究能否從中分析出含有這兩種頻率的信號

clear all %清除所有變量
N=256;dt=0.02;  %數據的個數和采樣間隔
n=0:N-1;t=n*dt; %序號序列和時間序列
x=sin(2*pi*t)+0.5*sin(2*pi*5*t);%合成的信號
m=floor(N/2)+1;%分解a,b的最大序列號，為分解的N/2個參數再加參數a0
%floor函數為向下取整
a=zeros(1,m);b=zeros(1,m);%產生a,b兩個為0的序列
for k=0:m-1
    for ii=0:N-1
        a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N);%求和
        b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N);%求和
    end
    %數組下標只能從1開始
    c(k+1)=sqrt(a(k+1).^2+b(k+1).^2);%k次諧波的振幅大小
end
subplot(2,1,1),plot(t,x);title('原始信號');xlabel('時間/s');ylim([-1.5 1.5])
%繪出時間域信號
subplot(2,1,2),plot((0:m-1)/(N*dt),c);
%繪出頻率域信號
title('Fourier變換');xlabel('頻率/Hz'),ylabel('振幅')

 可以看出Fourier變換識別出了1Hz和5Hz的波 與原信號不同是因為采樣點數較少

• 將m改為N，增加鏡像對稱的部分

clear all %清除所有變量
N=256;dt=0.02;  %數據的個數和采樣間隔
n=0:N-1;t=n*dt; %序號序列和時間序列
x=sin(2*pi*t)+0.5*sin(2*pi*5*t);%合成的信號
m=N;%鏡像對稱
%floor函數為向下取整
a=zeros(1,m);b=zeros(1,m);%產生a,b兩個為0的序列
for k=0:m-1
    for ii=0:N-1
        a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N);%求和
        b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N);%求和
    end
    %數組下標只能從1開始
    c(k+1)=sqrt(a(k+1).^2+b(k+1).^2);%k次諧波的振幅大小
end
subplot(2,1,1),plot(t,x);title('原始信號');xlabel('時間/s');ylim([-1.5 1.5])
%繪出時間域信號
subplot(2,1,2),plot((0:m-1)/(N*dt),c);
%繪出頻率域信號
title('Fourier變換');xlabel('頻率/Hz'),ylabel('振幅')

        

• m改為2N 頻譜增加一個周期

clear all %清除所有變量
N=256;dt=0.02;  %數據的個數和采樣間隔
n=0:N-1;t=n*dt; %序號序列和時間序列
x=sin(2*pi*t)+0.5*sin(2*pi*5*t);%合成的信號
m=2*N;%增加一個周期
%floor函數為向下取整
a=zeros(1,m);b=zeros(1,m);%產生a,b兩個為0的序列
for k=0:m-1
    for ii=0:N-1
        a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N);%求和
        b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N);%求和
    end
    %數組下標只能從1開始
    c(k+1)=sqrt(a(k+1).^2+b(k+1).^2);%k次諧波的振幅大小
end
subplot(2,1,1),plot(t,x);title('原始信號');xlabel('時間/s');ylim([-1.5 1.5])
%繪出時間域信號
subplot(2,1,2),plot((0:m-1)/(N*dt),c);
%繪出頻率域信號
title('Fourier變換');xlabel('頻率/Hz'),ylabel('振幅')

所以離散有限信號的頻譜為周期譜

• 根據分解得到的ak,bk恢復原來的信號

clear all %清除所有變量
N=256;dt=0.02;  %數據的個數和采樣間隔
n=0:N-1;t=n*dt; %序號序列和時間序列
x=sin(2*pi*t)+0.5*sin(2*pi*5*t);%合成的信號
m=floor(N/2)+1;%增加一個周期
%floor函數為向下取整
a=zeros(1,m);b=zeros(1,m);%產生a,b兩個為0的序列
for k=0:m-1
    for ii=0:N-1
        a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N);%求和
        b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N);%求和
    end
    %數組下標只能從1開始
    c(k+1)=sqrt(a(k+1).^2+b(k+1).^2);%k次諧波的振幅大小
end

if(mod(N,2)~=1)
    a(m)=a(m)/2;
end
%也就是N為偶數的時候 b(m)為0,a(m)減半
for ii=0:N-1
    xx(ii+1)=a(1)/2;
    for k=1:m-1
        xx(ii+1)=xx(ii+1)+a(k+1)*cos(2*pi*k*ii/N)+b(k+1)*sin(2*pi*k*ii/N);%求和
    end
end
subplot(2,1,1)
plot((0:N-1)*dt,x)%繪制原始信號
title('原始信號')
ylim([-2,2])
subplot(2,1,2)
plot((0:N-1)*dt,xx)%繪制合成信號
title('合成信號')
xlabel('時間/s')
ylim([-2,2])