%Samp3_2
clear all   %清除內(nèi)存所有變量
close all    %關閉所有打開的圖形窗口
N=200;dt=4/N; %數(shù)據(jù)點數(shù)和采樣間隔
for n=1:N    %得到方波信號
    if (n*dt>=2)
        x(n)=0.8;
    else
        x(n)=-0.8;
    end
end
figure(1)   %打開第一個圖形窗口，繪制原始信號及其Fourier分析系數(shù)
subplot(2,1,1),plot((1:N)*dt,x),hold on;
plot((1:N)*dt,zeros(1,N),'k'),xlabel('時間/s')
title('原始信號')
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('振幅') %繪制振幅譜
m=input('輸入諧波最大階數(shù)?');   %輸入最大k值
if(m>(floor(N/2)+1))   %如果最大k值大于Nyquist頻率對應的點數(shù)，則顯示出錯信息
    error('諧波最大階數(shù)必須小于Nyquist頻率對應的階數(shù).')
end
if(mod(N,2)~=1)a(nn)=a(nn)/2; end   %此時b(nn)為零，a(nn)減半
for ii=0:N-1    %合成信號
    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
figure(2)    %打開第二個圖形窗口，繪制合成信號圖
plot((1:N)*dt,xx,(0:N-1)*dt,x)    %繪制合成信號和方波圖形便于比較
hold on
plot((1:N)*dt,zeros(1,N),'k'),xlabel('時間/s')   %繪出橫軸
title('合成信號')