ultimatehomework.m

Matlab程序为清华大学胡广书版信号处理的上机习题答案

%------------------------------------------------------------------------
% ultimatehomework.m 最后大作业;
%------------------------------------------------------------------------
clear;

N=4096;
n=0:1:N-1;

%产生两个零均值、均匀分布、功率为0.02的白噪信号
sqer=0.02;
u1=rand(1,N);
u1=u1-mean(u1);
u1=u1*sqrt(12*sqer);       %设定u1的功率为0.02
u2=rand(1,N);
u2=u2-mean(u2);
u2=u2*sqrt(12*sqer);       %设定u2的功率为0.02

%FIR系统的5个级联的FIR子系统
b1=[1, 1.98,  0.9801];
b2=[1,-1.98,  0.9801];
b3=[1,-1.8418,0.9801];
b4=[1,-1.5,   0.9801];
b5=[1,-1.2727,0.9801];

%将信号u1和u2通过FIR滤波器
w1=filter(b1,1,u1);w1=filter(b2,1,w1);w1=filter(b3,1,w1);w1=filter(b4,1,w1);
v1=filter(b5,1,w1);clear w1;
w2=filter(b1,1,u2);w2=filter(b2,1,w2);w2=filter(b3,1,w2);w2=filter(b4,1,w2);
v2=filter(b5,1,w2);clear w2;
vn=v1+i*v2;

%四个复正弦参数
A1=15;A2=15;A3=7.5;A4=2;
fn1=0.14;fn2=0.15;fn3=0.24;fn4=-0.17;
zsin1=A1*exp(i*2*pi*fn1*n);
zsin2=A2*exp(i*2*pi*fn2*n);
zsin3=A3*exp(i*2*pi*fn3*n);
zsin4=A4*exp(i*2*pi*fn4*n);
xn=vn+zsin1+zsin2+zsin3+zsin4;

%求解信噪比
power_zsin1=dot(zsin1,zsin1)/N;
power_zsin2=dot(zsin2,zsin2)/N;
power_zsin3=dot(zsin3,zsin3)/N;
power_zsin4=dot(zsin4,zsin4)/N;
SNR1=10*log10(power_zsin1/(2*sqer));
SNR2=10*log10(power_zsin2/(2*sqer));
SNR3=10*log10(power_zsin3/(2*sqer));
SNR4=10*log10(power_zsin4/(2*sqer));

%作FIR系统的归一化幅频响应曲线
b=conv(b5,conv(b4,conv(b3,conv(b2,b1))));
[H_FIR,w]=freqz(b,1,N,'whole',1);
H_FIR=[H_FIR(N/2+1:N);H_FIR(1:N/2)];  %使在归一化频率-0.5~0.5上分布
w=[w(N/2+1:N)-1;w(1:N/2)];            %调整归一化频率的范围
Hr=abs(H_FIR);
Hr_max=max(Hr);                       
Hr_log=10*log10(Hr/Hr_max);
figure(1);
plot(w,Hr_log,'k');axis([min(w),0.5,-50,0]);grid on;
ylabel('|H[exp(jw)]| /dB');
xlabel('fs');

%作xn实部和虚部的波形图
xn_real=real(xn(1:256));
xn_imag=imag(xn(1:256));
figure(2);
subplot(211);
plot(n(1:256),xn_real,'k');axis([0,255,-40,40]);grid on;
ylabel('Re[x(n)]');
xlabel('n');
subplot(212);
plot(n(1:256),xn_imag,'k');axis([0,255,-40,40]);grid on;
ylabel('Im[x(n)]');
xlabel('n');

%求真实功率谱
power_real=2*sqer*(Hr.*Hr);
for wi=1:1:N
    if fn1>=w(wi) 
        wscr1=wi;
    end
    if fn2>=w(wi)
        wscr2=wi;
    end
    if fn3>=w(wi)
        wscr3=wi;
    end
    if fn4>=w(wi)
        wscr4=wi;
    end
end
power_real(wscr1)=power_real(wscr1)+power_zsin1;
power_real(wscr2)=power_real(wscr2)+power_zsin2;
power_real(wscr3)=power_real(wscr3)+power_zsin3;
power_real(wscr4)=power_real(wscr4)+power_zsin4;
power_real_max=max(power_real);                       
power_real_log=10*log10(power_real/power_real_max);

figure(3);
plot(w,power_real_log,'k');axis([min(w),0.5,-50,0]);grid on;
ylabel('Px(w) /dB  真实功率谱');
xlabel('fs');


%用周期图求功率谱
Nper=256;
wper=-0.5:1/N:0.5-1/N;
xn256=xn(1:Nper);%截取信号前256点
%xn256_amp=zeros(1,N);
%for wi=1:1:N
%    for Nper_i=1:1:Nper
%        xn256_amp(wi)=xn256_amp(wi)+...
%        xn256(Nper_i)*exp(-j*2*pi*wper(wi)*(Nper_i-1));
%    end
%end
    
%power_per=abs(xn256_amp).*abs(xn256_amp)/Nper;
%power_per_max=max(power_per);                       
%power_per_log=10*log10(power_per/power_per_max);
%figure(4);
%plot(wper,power_per_log,'k');axis([min(wper),0.5,-50,0]);grid on;
%ylabel('Pper(w) /dB  周期图功率谱');
%xlabel('fs');

%用Welch平均法估计信号的功率谱
power_welch=pwelch(xn256,hamming(64),32,N,'whole');        
power_welch_max=max(power_welch);
power_welch=power_welch/power_welch_max;
power_welch=10*log10(power_welch+0.000001);
figure(5);
plot(wper,fftshift(power_welch),'k');axis([min(wper),0.5,-50,0]);grid on;
ylabel('Pper(w) /dB  Welch平均法功率谱');
xlabel('fs');

%用AR模型估计信号的功率谱自相关法和Burg算法
power_yulear=pyulear(xn256,10,N);
power_yulear_max=max(power_yulear);
power_yulear=power_yulear/power_yulear_max;
power_yulear=10*log10(power_yulear+0.000001);
figure(6);
plot(wper,fftshift(power_yulear),'k');axis([min(wper),0.5,-50,0]);grid on;
ylabel('Par(w) /dB  AR模型自相关功率谱');
xlabel('fs');

power_burg=pburg(xn256,10,N);
power_burg_max=abs(max(power_burg));
power_burg=abs(power_burg)/power_burg_max;
power_burg=10*log10(power_burg+0.000001);
figure(7);
plot(wper,fftshift(power_burg),'k');axis([min(wper),0.5,-50,0]);grid on;
ylabel('Par(w) /dB  AR模型Burg功率谱');
xlabel('fs');

%估计自相关法和Burg算法的阶次
%orderM=40;
%eni_yule=zeros(1,orderM);
%for mi=1:1:orderM
%    [are,eni_yule(mi)]=aryule(xn256,mi);
%end
%eni_yule=abs(eni_yule);
%eni_burg=zeros(1,orderM);
%for mi=1:1:orderM
%    [are,eni_burg(mi)]=arburg(xn256,mi);
%end
%eni_burg=abs(eni_burg);

%km=1:1:orderM;
%FPE_yule=eni_yule.*(Nper+km+1)./(Nper-km-1);
%FPE_burg=eni_burg.*(Nper+km+1)./(Nper-km-1);
%figure(8);
%subplot(221);plot(km,FPE_yule,'k');
%axis([1,orderM,0,100]);grid on;
%ylabel('FPE　  AR模型Yule阶次');xlabel('p  阶次');
%subplot(222);plot(km,FPE_burg,'k');
%axis([1,orderM,0,100]);grid on;
%ylabel('FPE　  AR模型Burg阶次');xlabel('p  阶次');

% 设计切比雪夫最佳一致逼近滤波器；
order_flt=56;
fps=[0 0.20 0.23 0.5];fps=2*fps;
Amp=[1   1   0    0 ];
weigh=[2.72 1];	
b_chv_flt=remez(order_flt,fps,Amp,weigh);

[h_chv_flt,w_chv_flt]=freqz(b_chv_flt,1,256,'whole',1);
hr_chv_flt=abs(h_chv_flt);
hr_chv_flt=hr_chv_flt/max(hr_chv_flt);%幅度归一化
hp_chv_flt=unwrap(angle(h_chv_flt));
hrlg_chv_flt=20*log10(hr_chv_flt);

figure(9)
stem(b_chv_flt,'k.');axis([0,60,-0.1,0.501]);grid;ylabel('h(n)');xlabel('n');
figure(10)
plot(w_chv_flt-0.5,fftshift(hrlg_chv_flt),'k');axis([-0.5,0.5,-50,10]);grid on;
ylabel('|H[exp(jw)]| /dB');xlabel('fs');
figure(11)
plot(w_chv_flt(1:128),hp_chv_flt(1:128),'k');grid on;
ylabel('hpase');xlabel('fs');

yn=filter(b_chv_flt,1,xn);
yn256=yn(1:256);

%作yn实部和虚部的波形图
yn_real=real(yn256(1:256));
yn_imag=imag(yn256(1:256));
figure(12);
subplot(211);
plot(n(1:256),yn_real,'k');axis([0,256,-40,40]);grid on;
ylabel('Re[y(n)]');
xlabel('n');
subplot(212);
plot(n(1:256),yn_imag,'k');axis([0,255,-40,40]);grid on;
ylabel('Im[y(n)]');
xlabel('n');

%用周期图求yn功率谱
Nper=256;
wper=-0.5:1/N:0.5-1/N;
yn256=yn(1:Nper);%截取信号前256点
%yn256_amp=zeros(1,N);
%for wi=1:1:N
%    for Nper_i=1:1:Nper
%        yn256_amp(wi)=yn256_amp(wi)+...
%        yn256(Nper_i)*exp(-j*2*pi*wper(wi)*(Nper_i-1));
%    end
%end
    
%power_per=abs(yn256_amp).*abs(yn256_amp)/Nper;
%power_per_max=max(power_per);                       
%power_per_log=10*log10(power_per/power_per_max);
%figure(13);
%plot(wper,power_per_log,'k');axis([min(wper),0.5,-50,0]);grid on;
%ylabel('Pper(w) /dB  y(n)周期图功率谱');
%xlabel('fs');

%用Welch平均法估计信号的功率谱
power_welch=pwelch(yn256,hamming(64),32,N,'whole');        
power_welch_max=max(power_welch);
power_welch=power_welch/power_welch_max;
power_welch=10*log10(power_welch+0.000001);
figure(14);
plot(wper,fftshift(power_welch),'k');axis([min(wper),0.5,-50,0]);grid on;
ylabel('Pper(w) /dB  Welch平均法y(n)功率谱');
xlabel('fs');

%用AR模型估计信号的功率谱自相关法和Burg算法
power_yulear=pyulear(yn256,25,N);
power_yulear_max=max(power_yulear);
power_yulear=power_yulear/power_yulear_max;
power_yulear=10*log10(power_yulear+0.000001);
figure(17);
plot(wper,fftshift(power_yulear),'k');axis([min(wper),0.5,-50,0]);grid on;
ylabel('Par(w) /dB  AR模型自相关功率谱');
xlabel('fs');

power_burg=pburg(yn256,25,N);
power_burg_max=abs(max(power_burg));
power_burg=abs(power_burg)/power_burg_max;
power_burg=10*log10(power_burg+0.000001);
figure(18);
plot(wper,fftshift(power_burg),'k');axis([min(wper),0.5,-50,0]);grid on;
ylabel('Par(w) /dB  AR模型Burg功率谱');
xlabel('fs');

%估计自相关法和Burg算法的阶次
orderM=40;
eni_yule=zeros(1,orderM);
for mi=1:1:orderM
    [are,eni_yule(mi)]=aryule(yn256,mi);
end
eni_yule=abs(eni_yule);
eni_burg=zeros(1,orderM);
for mi=1:1:orderM
    [are,eni_burg(mi)]=arburg(yn256,mi);
end
eni_burg=abs(eni_burg);

km=1:1:orderM;
FPE_yule=eni_yule.*(Nper+km+1)./(Nper-km-1);
FPE_burg=eni_burg.*(Nper+km+1)./(Nper-km-1);
figure(19);
subplot(221);plot(km,FPE_yule,'k');
axis([1,orderM,0,15]);grid on;
ylabel('FPE　  AR模型Yule阶次');xlabel('p  阶次');
subplot(222);plot(km,FPE_burg,'k');
axis([1,orderM,0,15]);grid on;
ylabel('FPE　  AR模型Burg阶次');xlabel('p  阶次');

%Pisarenko谐波分解法
rowN=5;colN=5;
acr_xn=xcorr(xn256,rowN,'biased');
acr_xn=acr_xn(rowN+1:2*rowN+1);
Rmatrix=zeros(rowN,colN);        %创建自相关矩阵　
for rowi=1:1:rowN
    for coli=1:1:rowi
        Rmatrix(rowi,coli)=acr_xn(rowi-coli+1);
    end
    for coli=rowi+1:1:colN
        Rmatrix(rowi,coli)=acr_xn(coli-rowi+1)';
    end
end
[eivector,eivalue]=eig(Rmatrix);
vm1=eivector(:,1:1);
rt=roots(vm1);
omg=angle(rt)/2/pi;
Wmatrix=zeros(rowN-1,colN-1);
for rowi=1:1:rowN-1
    for coli=1:1:colN-1
        Wmatrix(rowi,coli)=exp(j*rowi*omg(coli));
    end
end
Am=(Wmatrix)^(-1)*acr_xn(1:rowN-1)'%振幅的平方，即功率的计算

[power_music,fw]=pmusic(xn256,10,4096,1);
power_music_max=max(power_music);
power_music=power_music/power_music_max;
power_music=10*log10(power_music+0.000001);
figure(20);
plot(fw-0.5,fftshift(power_music),'k');axis([-0.5,0.5,-50,0]);grid on;
ylabel('Par(w) /dB  非参数模型功率谱');
xlabel('fs');




%fn=fopen('result35gs7.dat','w');
%for i=1:1:M
%    fprintf(fn,'%12.2f    %12.8f\n',n2(i),logP(i));
%end
%fclose(fn);


%[H1,w1]=freqz(b1,1,N,'whole',1);
%[H2,w]=freqz(b2,1,N,'whole',1);
%[H3,w]=freqz(b3,1,N,'whole',1);
%[H4,w]=freqz(b4,1,N,'whole',1);
%[H5,w]=freqz(b5,1,N,'whole',1);
%subplot(211);