• 將振幅為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])