现代数字信号matlab处理仿真题.docVIP

PAGE 2 3.17 （1）相关函数 仿真代码： A1=getAk(SNR1); A2=getAk(SNR2); A3=getAk(SNR3); %求得信号的幅度； noise=randn(1,N) + j*randn(1,N); % 构建高斯白噪声； s1=getSk(A1,f1,N); s2=getSk(A2,f2,N); s3=getSk(A3,f3,N);; %产生3个复正弦信号 vn=s1+s2+s3+noise; vk=fft(vn,2*N); %对v(n)补N个零，然后做2N点FFT swk=((abs(vk)).^2)/N; %计算功率谱估计S(ω) r0=ifft(swk); %对S(k)做ifft得到 r=[r0(N+2 : 2*N) , r0(1 : N)]; %根据教程3.1.8式可得 r1=xcorr(vn, N-1,biased); %直接计算自相关函数 %%%%%%%%%%%%%%%%%%%%%%%%%取序列实部，虚部%%%%%% real_r=real(r); imag_r=imag(r); real_r1=real(r1); imag_r1 = imag(r1); subplot(2,2,1); stem(real_r); xlabel(基于FFT的自相关函数的实部); ylabel(实部); subplot(2,2,2); stem(imag_r); xlabel(基于FFT的自相关函数的虚部); ylabel(虚部); subplot(2,2,3); stem(real_r1); ylabel(实部); xlabel(估计的自相关函数的实部); subplot(2,2,4); stem(imag_r1); xlabel(估计的自相关函数的虚实部); ylabel(虚部); function AK=getAk(SNR) %求得幅度 %%%%%%%%%%%%%%%%%%%由SNR=10log(A^2/2*σ^2) %%%%%%% AK=((10^(SNR/10))*2)^0.5; function Sk=getSk(Ak,fk,N) Sk=Ak * exp(j * 2 * pi * fk *(0:N-1)); 仿真波形： （2）BT法和周期法估计 仿真程序： clear all; clc; %设定N值可以改变抽样信号的点数，设定M值可以设定加窗的大小，设定N3可以补零，确定实际求fft的点数。 N=256; %样本观察长度 M=64; %加窗长度 f1=0.15; f2=0.17; f3=0.26; %定义3个归一化正弦频率 SNR1=30; SNR2=30; SNR3=27; %定义三个正弦信号信噪比 A1=getAk(SNR1); A2=getAk(SNR2); A3=getAk(SNR3); %求得信号的幅度； noise=randn(1,N) + j *randn(1,N); % 构建高斯白噪声； s1=getSk(A1,f1,N); s2=getSk(A2,f2,N); s3=getSk(A3,f3,N); %产生3个复正弦信号 un=s1+s2+s3+noise; %%%%%%%%%%%%下面采用周期法估计的频谱%%%%%%%%%%%%% N2=2*N; u2=zeros(1,N2); u2(1:N)=un;%对信号序列进行补零 Uw=fft(u2);%求DFT变换 Sw1=zeros(1,N2); for w=1:N2 Sw1(w)=((abs(Uw(w)))^2)/N;%计算功率谱 end r0=zeros(1,N2); r0=ifft(Sw1); r=zeros(1,N2); r(1:N)=r0(N+1:N2); r(N+1:N2)=r0(1:N); M=64;%加窗处理 for n=1:2*M r1(n)=r((N2-2*M)/2+n); end%加窗之后的序列为r1 N3=256;%实际进行fft变换的点数。 r2=zeros(1,N3); r2(1:2*M)=r1;%将加窗之后的r1进行了补零得到的是r2. Sw2=fft(r2); Sw2=log10(abs(Sw2));%取模取对数 temp=zeros(1,N3/2); temp=Sw2(1:N3/2); Sw2(1:N3/2)=Sw2((N3/2+1):N3); Sw2((N3/2+1):N3)=temp;%将0-2*pi的序列折换成-pi到pi的序列。 d=1/(N3-1); x=zeros(1,N3); for i=1:N3 x(i)=-0.5+(i-1)*d; end Sw1=log10(abs(Sw1));%取模取对数 temp=z