%  ex_12_11_1.m :  習(xí)題 12.11


clear all;
close all;

%題目給出的參數(shù)-------------------------------------------------------
N=256;
P=0.12;     
h=[1 -0.1 0.09 0.648];        %FIR濾波器系數(shù)
SNR=[10 50 50];   F=[0.1 0.25 0.26];          %三個信噪比和歸一化頻率

%步驟一：產(chǎn)生白噪聲-----------------------------------------------
a=sqrt(P);
u=rand(1,20*N); 
u=u-mean(u);u=u*a;

%步驟二：濾波處理-----------------------------------------------------
v1=filter(h,1,u);
v=v1(1000:1000+N-1);

%步驟三：產(chǎn)生正弦序列----------------------------------------------------
A=sqrt(2*P*10.^(SNR/10))    %根據(jù)信噪比確定三個正弦函數(shù)的幅度
x1=A(1)*sin(2*pi*F(1)*[0:N-1]);
x2=A(2)*sin(2*pi*F(2)*[0:N-1]);
x3=A(3)*sin(2*pi*F(3)*[0:N-1]);

%步驟四：疊加信號，得到信號x---------------------------------------------
x=v+x1+x2+x3;
save testdat x;

%步驟五：求出信號x的真實功率譜---------------------------------------------
M=4096;H=freqz(h,1,M,1);H=abs(H);H=H.*H;H=P*H;
f=0:1/M:0.5-1/M;
power_y1=pi*A(1)*A(1)/2;
power_y2=pi*A(2)*A(2)/2;
power_y3=pi*A(3)*A(3)/2;
f1=0.1;f2=0.25;f3=0.26;
f1_n=round(f1*M);
f2_n=round(f2*M);
f3_n=round(f3*M);

H(f1_n)=H(f1_n)+power_y1;
H(f2_n)=H(f2_n)+power_y2;
H(f3_n)=H(f3_n)+power_y3;
PSD=10*log10(H/max(H)+eps);
figure('color','w');plot(f,PSD(1:2048));grid on;