%Samp3_2
clear all   %清除內(nèi)存所有變量
close all    %關(guān)閉所有打開的圖形窗口
N=200;dt=4/N; %數(shù)據(jù)點(diǎn)數(shù)和采樣間隔
for n=1:N    %得到方波信號(hào)
    if (n*dt>=2)
        x(n)=0.8;
    else
        x(n)=-0.8;
    end
end
figure(1)   %打開第一個(gè)圖形窗口，繪制原始信號(hào)及其Fourier分析系數(shù)
subplot(2,1,1),plot((1:N)*dt,x),hold on;
plot((1:N)*dt,zeros(1,N),'k'),xlabel('時(shí)間/s')
title('原始信號(hào)')
a=zeros(1,N);b=zeros(1,N);
nn=floor(N/2)+1;
for k=0:nn-1
    a(k+1)=0;
    b(k+1)=0;
    for ii=0:N-1
        a(k+1)=a(k+1)+2/N*x(ii+1)*cos(2*pi*k*ii/N);  %求解Fourier系數(shù)
        b(k+1)=b(k+1)+2/N*x(ii+1)*sin(2*pi*k*ii/N);
    end
    c(k+1)=sqrt(a(k+1).^2+b(k+1).^2);
end
subplot(2,1,2),plot((0:nn-1)/(N*dt),c);
title('Fourier變換')
xlabel('頻率/Hz'),ylabel('振幅') %繪制振幅譜
figure(2)    %打開第二個(gè)圖形窗口，繪制合成信號(hào)圖
if(mod(N,2)~=1)a(nn)=a(nn)/2; end  %此時(shí)b(nn)為零，a(nn)減半
m=1;
for ii=0:N-1    %合成信號(hào)
    xx(ii+1)=a(1)/2;
    for k=1:m
        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,2,1)
plot((1:N)*dt,xx,(0:N-1)*dt,x)    %繪制合成信號(hào)和方波圖形便于比較
hold on
plot((1:N)*dt,zeros(1,N),'k'),xlabel('時(shí)間/s')   %繪出橫軸
title('m=1的合成信號(hào)')

m=3;
for ii=0:N-1    %合成信號(hào)
    xx(ii+1)=a(1)/2;
    for k=1:m
        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,2,2)
plot((1:N)*dt,xx,(0:N-1)*dt,x)    %繪制合成信號(hào)和方波圖形便于比較
hold on
plot((1:N)*dt,zeros(1,N),'k'),xlabel('時(shí)間/s')   %繪出橫軸
title('m=3的合成信號(hào)')

m=5;
for ii=0:N-1    %合成信號(hào)
    xx(ii+1)=a(1)/2;
    for k=1:m
        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,2,3)
plot((1:N)*dt,xx,(0:N-1)*dt,x)    %繪制合成信號(hào)和方波圖形便于比較
hold on
plot((1:N)*dt,zeros(1,N),'k'),xlabel('時(shí)間/s')   %繪出橫軸
title('m=5的合成信號(hào)')

m=7;
for ii=0:N-1    %合成信號(hào)
    xx(ii+1)=a(1)/2;
    for k=1:m
        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,2,4)
plot((1:N)*dt,xx,(0:N-1)*dt,x)    %繪制合成信號(hào)和方波圖形便于比較
hold on
plot((1:N)*dt,zeros(1,N),'k'),xlabel('時(shí)間/s')   %繪出橫軸
title('m=7的合成信號(hào)')