expfft.m

抽样定理的实验体会

function fft(xx, yy, zz, aa)
if zz==0
    p=xx;
    q=yy;
    for n=1:16
        x(n)=exp(-(n-1-p)^2/q);
    end
    subplot(2,1,1);
    stem((0:15),x(1:16));
    grid;
    xlabel('N');
    ylabel('幅度');
    title('Xn');
    %f
    subplot(2,1,2);
    h(1:16)=fft(x(1:16));
    stem((0:15),abs(h(1:16)));
    xlabel('K');
    ylabel('幅度');
    title('|Xk|');
    grid;   
elseif zz==1
    a=xx; 
    f=yy;
    for n=1:16
        x(n)=exp(-a*(n-1))*sin(2*pi*f*(n-1));
    end
    subplot(2,1,1);
    stem((0:15),x(1:16));
    grid;
    xlabel('N');
    ylabel('幅度');
    title('Xn');
    %f
    subplot(2,1,2);
    h(1:16)=fft(x(1:16));
    stem((0:15),abs(h(1:16)));
    grid; 
    xlabel('K');
    ylabel('幅度');
    title('|Xk|');
elseif zz==2
    for n=0:3
        xc(n+1)=n;
    end
    for n=4:7
        xc(n+1)=8-n;
    end
    for n=9:16
        xc(n)=0;
    end
    for n=0:3
       xd(n+1)=4-n;
    end
    for n=4:7
       xd(n+1)=n-4;
    end
    for n=9:16
       xd(n)=0;
    end
    if aa==0
        subplot(211),stem((0:7),xc(1:8));
        grid;
        title('8点正三角波波形');
        subplot(212),stem((0:7),xd(1:8));
        grid;
        title('8点反三角波波形');
    elseif aa==1
        hc(1:8)=fft(xc(1:8));
        subplot(211),stem((0:7),abs(hc(1:8)));
        grid;
        title('8点正三角波的幅频特性');
        hd(1:8)=fft(xd(1:8));
        subplot(212),stem((0:7),abs(hd(1:8)));
        grid;
        title('8点反三角波的幅频特性');
    elseif aa==2
        subplot(211),stem((0:15),xc(1:16));
        title('8点正三角波补0到16点后的时间波形');
        grid;
        subplot(212),stem((0:15),xd(1:16));
        title('8点反三角波补0到16点后的时间波形');
        grid;
    elseif aa==3
        hc(1:16)=fft(xc(1:16));
        subplot(211),stem((0:15),abs(hc(1:16)));
        title('8点正三角波补0到16点后的幅频特性');
        grid;
        hd(1:16)=fft(xd(1:16));
        subplot(212),stem((0:15),abs(hd(1:16)));
        title('8点反三角波补0到16点后的幅频特性');
        grid;
    end
elseif zz==3
    N=16;
    deltf=1/16;
    for n=1:N
        x(n)=sin(2*pi*0.125*n)+cos(2*pi*(0.125+deltf)*n);
    end
    deltf=1/64;
    for n=1:N
        x1(n)=sin(2*pi*0.125*n)+cos(2*pi*(0.125+deltf)*n);
    end
    if aa==0
        subplot(211),stem((0:N-1),x(1:N));
        grid;
        xlabel('N');
        ylabel('幅度');
        title('Xn(N=16,deltf=1/16)');
        subplot(212),stem((0:N-1),x1(1:N));
        grid;
        xlabel('N');
        ylabel('幅度');
        title('Xn(N=16,deltf=1/64)');
    elseif aa==1
        h(1:N)=fft(x(1:N));
        subplot(211),stem((0:N-1),abs(h(1:N)));
        xlabel('K');
        ylabel('幅度');
        title('|Xk|(N=16,deltf=1/16)');
        grid;
        h1(1:N)=fft(x1(1:N));
        subplot(212),stem((0:N-1),abs(h1(1:N)));
        xlabel('K');
        ylabel('幅度');
        title('|Xk|(N=16,deltf=1/64)');
        grid;        
    elseif aa==2
        N=128;
        deltf=1/16;
        for n=1:N
            x(n)=sin(2*pi*0.125*n)+cos(2*pi*(0.125+deltf)*n);
        end
        subplot(211),stem((0:N-1),x(1:N));
        grid;
        xlabel('N');
        ylabel('幅度');
        title('Xn(N=128,deltf=1/16)');
        deltf=1/64;
        for n=1:N
            x1(n)=sin(2*pi*0.125*n)+cos(2*pi*(0.125+deltf)*n);
        end
        subplot(212),stem((0:N-1),x1(1:N));
        grid;
        xlabel('N');
        ylabel('幅度');
        title('Xn(N=128,deltf=1/64)');
    elseif aa==3
        h(1:N)=fft(x(1:N));
        subplot(211),stem((0:N-1),abs(h(1:N)));
        xlabel('K');
        ylabel('幅度');
        title('|Xk|(N=128,deltf=1/16)');
        grid;
        h1(1:N)=fft(x1(1:N));
        subplot(212),stem((0:N-1),abs(h1(1:N)));
        xlabel('K');
        ylabel('幅度');
        title('|Xk|(N=128,deltf=1/64)');
        grid;
    end
elseif zz==4
    p=8;
    q=2;
    for n=1:16
       xa(n)=exp(-(n-1-p)^2/q);
    end
    for n=17:32
      xa(n)=0;
    end
    a=0.1; 
    f=0.0625;
    for n=1:16
        xb(n)=exp(-a*(n-1))*sin(2*pi*f*(n-1));
    end
    for n=17:32
       xb(n)=0;
    end
    if aa==0
        subplot(211),stem((0:15),xa(1:16));
        grid;
        xlabel('N');
        ylabel('幅度');
        title('xa(n)=exp(-(n-8)^2/2)');
        subplot(212),stem((0:15),xb(1:16));
        grid;
        xlabel('N');
        ylabel('幅度');
        title('xb(n)=exp(-a*n)*sin(2*pi*f*n)');
    elseif aa==1
        ha(1:16)=fft(xa(1:16));
        hb(1:16)=fft(xb(1:16));
        h1(1:16)=ha(1:16).*hb(1:16);
        f1(1:16)=ifft(h1(1:16));
        ha(1:32)=fft(xa(1:32));
        hb(1:32)=fft(xb(1:32));
        h2(1:32)=ha(1:32).*hb(1:32);
        f2(1:32)=ifft(h2(1:32));
        subplot(211),stem((0:15),real(f1(1:16)));
        grid;
        xlabel('N');
        ylabel('幅度');
        title('圆周卷积（N=16)：f1(n)=xa(n)*xb(n)');
        subplot(212),stem((0:31),real(f2(1:32)));
        grid;
        xlabel('N');
        ylabel('幅度');
        title('线性卷积（N=32)：f2(n)=xa(n)*xb(n)          注：在n=10~15点与圆周卷积相等');
    end
elseif zz==5
    N=512;
    for n=1:4
        xc(n)=n;
    end
    for n=5:8
        xc(n)=9-n;
    end
    for n=9:16
        xc(n)=0;
    end
    xe=randn(1,N);
    %exp6_0 plot    
    %exp6_1 plot
    
    
        %重叠相加法
        L=64;
        M=ceil(N/L);%进一法取整数
        f1=zeros(1,M*L+16-1);
        xe=[xe,zeros(1,M*L-N)];%如果xe不够L的整数倍长度，则补零
        L2=L*2;%分段卷积时数据长度要扩展，为使用FFT算法，应扩展到2的整数次方
        x1=[xc,zeros(1,L2-16)];%数据长度扩展部分要补0
        X1=fft(x1);
        for i=0:M-1
            x=xe(1+L*i:L+L*i);
            x2=[x,zeros(1,L2-L)];%数据长度扩展部分要补0
            X2=fft(x2);
            f0=ifft(X1.*X2);%用FFT方法计算xe(n)中第i段序列与xc(n)的线性卷积
            f=real(f0(1:L+16-1));%根据公式推导，线性卷积只有实部，有效长度为L+16-1
            f1(1+L*i:L+16-1+L*i)=f1(1+L*i:L+16-1+L*i)+f;%将分段线性卷积值加到对应的段中
        end
        %重叠保留法
        L=64;
        M=ceil(N/L);%进一法取整数
        f2=zeros(1,M*L+16-1);
        xe=[xe,zeros(1,M*L-N)];%如果xe不够L的整数倍长度，则补零
        L2=L*2;%分段卷积时数据长度要扩展，为使用FFT算法，应扩展到2的整数次方
        x1=[xc,zeros(1,L2-16)];%数据长度扩展部分要补0
        X1=fft(x1);
        x2=[zeros(1,L2-L),xe(1:L)];%在xe(n)序列前补0
        X2=fft(x2);
        f0=ifft(X1.*X2);%用FFT方法计算xe(n)中第1段序列与xc(n)的线性卷积
        f2(1:L)=real(f0(L2-L+1:L2));%取后面的L个数据，并取实部
        for i=0:M-2
            x2=xe(1+L*i:L2+L*i);
            X2=fft(x2);
            f0=ifft(X1.*X2);%用FFT方法计算xe(n)中第i段序列与xc(n)的线性卷积
            f=real(f0(L2-L+1:L2));%取后面的L个数据，并取实部
            f2(1+L*(i+1):L+L*(i+1))=f2(1+L*(i+1):L+L*(i+1))+f;%循环外已有一段，从第2段起保存
        end
        x2=[xe(1+L*(M-1):L*M),zeros(1,L2-L)];%在xe(n)后补L2-L个0
        X2=fft(x2);
        f0=ifft(X1.*X2);%用FFT方法计算xe(n)中最后一段序列与xc(n)的线性卷积
        f=real(f0(L2-L+1:L2-L+16-1));%最后一段只需取xc(n)数据长度并可少1个数据
        f2(1+L*M:16-1+L*M)=f2(1+L*M:16-1+L*M)+f;%
        %直接用FFT作卷积
        x1=[xc,zeros(1,N-1)];%总长度为N+16-1，本身为16点，补N-1个0
        X1=fft(x1);
        x2=[xe,zeros(1,16-1)];%总长度为N+16-1，本身为N点，补16-1个0
        X2=fft(x2);
        F3=X1.*X2;
        f0=ifft(F3);%用FFT方法计算xe(n)中第i段序列与xc(n)的线性卷积
        f3(1:N+16-1)=real(f0(1:N+16-1));%取线性卷积有效长度，并取实部 
        %exp6_2 plot
        
     %exp6_3 plot
       
    if aa==0
        subplot(211),stem((0:15),xc(1:16));
        title('8点正三角波补0到16点后的时间波形');
        grid;
        subplot(212),plot((0:N-1),xe(1:N));
        grid;
        title('512点标准正态分布随机序列波形');
    elseif aa==1
        hc=fft(xc);
        subplot(211),stem((0:15),abs(hc(1:16)));
        title('8点正三角波补0到16点后的幅频特性');
        grid;
        he=fft(xe);
        subplot(212),plot((0:N-1),abs(he(1:N)));
        title('512点标准正态分布随机序列的幅频特性');
        grid;
    elseif aa==2    
        subplot(311),plot(f3);
        title('直接FFT计算法线性卷积结果');
        grid;
        subplot(312),plot(f1);
        title('重叠相加法线性卷积结果');
        grid;
        subplot(313),plot(f2);
        title('重叠保留法线性卷积结果');
        grid;
    elseif aa==3   
        F1=fft(f1);
        F2=fft(f2);
        subplot(311),plot(abs(F3));
        title('直接FFT计算法线性卷积结果');
        grid;
        subplot(312),plot(abs(F1));
        title('重叠相加法线性卷积结果');
        grid;
        subplot(313),plot(abs(F2));
        title('重叠保留法线性卷积结果');
        grid;
    end
    
end
clear;