设信号x(t)=cos(2π×50t)+2×cos(2π×400t)，试将它的两个频率分量分离，并绘制它们的时域波形及频谱图

以下程序无需赋值，直接运行即可：

（已验证可以运行）

function [yl,yh]=shiyan49

fs=1600; %采样频率

Tt=0.02; %信号周期

T0=4*Tt; %记录长度

[xn,wk,N]=shiyan40(fs,T0);

M=length(wk);

if M==2 

rp=1;rs=80;

    f1=wk(1)*fs/N;

    f2=wk(2)*fs/N;

    f0=f2-f1;

    fpl=f1+f0/20; fstl=f2-f0/10; %模拟低通滤波器的特征参数

    [bzl,azl]=shiyan42(fpl,fstl,rp,rs,fs); 

    fph=f2-f0/2;fsth=f1+f0/20; %模拟高通滤波器的特征参数

    [bzh,azh]=shiyan43(fph,fsth,rp,rs,fs); 

end

ynl=filter(bzl,azl,xn); %序列xn通过数字低通滤波器，输出为ynl

ynh=filter(bzh,azh,xn); %序列xn通过数字高通滤波器，输出为ynh

knl=abs(fft(ynl)); % ynl的频谱knl

kl=knl/max(knl); % ynl的幅度归一化频谱kl 

knh=abs(fft(ynh)); kh=knh/max(knh);

T=1/fs;

figure(8)

t=0:T:(T0-T);

w=0:2*pi/N:(2*pi-2*pi/N);

subplot(2,2,1); plot(t,ynl);

subplot(2,2,2); stem(w,kl);

subplot(2,2,3); plot(t,ynh);

subplot(2,2,4); stem(w,kh);

%去掉滤波输出序列的头一个周期

yl=ynl(N/4+1:N);

yh=ynh(N/4+1:N);

knl=abs(fft(yl)); kl=knl/max(knl);

knh=abs(fft(yh)); kh=knh/max(knh);

N=length(kl);

t=Tt:T:(T0-T);

w=0:2*pi/N:(2*pi-2*pi/N);

figure(9)

subplot(2,2,1)

plot(t,yl);

xlabel('秒'); title('低频分量时域波形')

subplot(2,2,2); 

stem(w,kl)

xlabel('数字频率/弧度 ');title('低频分量归一化频谱图')

subplot(2,2,3)

plot(t,yh);

xlabel('秒'); title('高频分量时域波形')

subplot(2,2,4)

stem(w,kh)

xlabel('数字频率/弧度 '); title('高频分量归一化频谱图')

 

function [xn,wk,N]=shiyan40(fs,T0)

T=1/fs;

t=0:T:(T0-T);

xn=cos(100*pi*t)+2*cos(800*pi*t);

xk=fft(xn);

N=length(xn);

i=1;

wk=0;

for m=1:1:(N+1)/2

   if (abs(xk(m))>0.0001)

       wk(i)=m-1;

       i=i+1;

    end

end

n=0:1:N-1;

figure(1)

subplot(2,1,1)

plot(t,xn);

subplot(2,1,2)

stem(n,abs(xk),'r')

 

function [Ba,Aa,Wn]=shiyan41(mp,ms,rp,rs) 

%生成归一化频率的模拟低通滤波器

[N,Wn]=buttord(mp,ms,rp,rs,'s');%mp/ms:通带/阻带截止频率(弧度%/秒)

[z,p,k]=buttap(N);

[Ba,Aa]=zp2tf(z,p,k);

 

function [bz,az]=shiyan42(fp,fst,rp,rs,fs) 

%数字低通滤波器的生成

W0=[0,fp,fst,fs/2]; rr=[rp,rp,rs,rs];%设计指标

mp=fp*2*pi;

ms=fst*2*pi; %mp/ms:通带/阻带截止频率（弧度/秒）

[Ba,Aa,Wn]=shiyan41(mp,ms,rp,rs);

[b,a]=lp2lp(Ba,Aa,Wn);%将归一化频率的低通滤波器转换成截止频

%率为Wn的数字低通滤波器

[bz,az]=bilinear(b,a,fs);

[H,W]=freqz(bz,az);   

figure(2);

plot(W*fs/(2*pi),20*log10(abs(H)));

hold on

plot(W0,-rr,'r')

xlabel('频率/Hz')

ylabel('幅频特性/dB')

grid

 

function [bz,az]=shiyan43(fp,fst,rp,rs,fs) 

%数字高通滤波器的生成

W0=[0,fst,fp,fs/2]; rr=[rs,rs,rp,rp]; %设计指标

wp=fp*2*pi/fs; ws=fst*2*pi/fs; %模拟指标数字化

    C=2*fs; %频率预畸

    hp=C*tan(wp/2); hs=C*tan(ws/2);

    mp=hp; ms=hp^2/hs; %模拟高通指标转化为模拟低通指标

[Ba,Aa,Wn]=shiyan41(mp,ms,rp,rs);

[b,a]=lp2hp(Ba,Aa,mp); %将归一化频率的低通滤波器转换成截止

%频率为Wn的高通滤波器

[bz,az]=bilinear(b,a,fs);

[H,W]=freqz(bz,az);

figure(3)

plot(W*fs/(2*pi),20*log10(abs(H)));

hold on

plot(W0,-rr,'r')

xlabel('频率/Hz')

ylabel('幅频特性/dB')

grid

 

 

运行结果截图：