function st=mstg
N=1600
Fs=10000;T=1/Fs;Tp=N*T;
t=0:T:(N-1)*T;k=0:N-1;f=k/Tp;
fc1=Fs/10;
fm1=fc1/10;
fc2=Fs/20;
fm2=fc2/10;
fc3=Fs/40;
fm3=fc3/10;
xt1=cos(2*pi*fm1*t).*cos(2*pi*fc1*t);
xt2=cos(2*pi*fm2*t).*cos(2*pi*fc2*t);
xt3=cos(2*pi*fm3*t).*cos(2*pi*fc3*t);
st=xt1+xt2+xt3;
fxt=fft(st,N);
subplot(3,1,1)
plot(t,st);grid;xlabel('t/s');ylabel('s(t)');
axis([0,Tp/8,min(st),max(st)]);title('(a)s(t)的波形')
subplot(3,1,2)
stem(f,abs(fxt)/max(abs(fxt)),'.');grid ;title('(b) s(t)的频谱')
axis([0,Fs/5,0,1.2]);
xlabel('f/Hz'),ylabel('幅度')
数字滤波器的设计
clear all;close all;
Fs=10000; T=1/Fs;
st=mstg;
fp=290;fs=440
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;
[N,wp]=ellipord(wp,ws,rp,rs);
[B,A]=ellip(N,rp,rs,wp);
y1t=filter(B,A,st);
figure(2);subplot(2,1,1);
[H,W]=freqz(B,A,1000);
m=abs(H);
plot(W/pi,20*log10(m));grid on;
xlabel('\\omega/\\pi');ylabel('-A(dB)')
axis([0,1,-110,5]);title('损耗函数曲线');
figure(2);subplot(2,1,2);
n=0:length(y1t)-1;t=n*T;
plot(t,y1t);grid on;
xlabel('t/s');ylabel('yn');
axis([0,t(end),min(y1t),1.2*max(y1t)]);title('输出波形');
clear all;close all;
Fs=10000; T=1/Fs;
st=mstg;
fp1=430;fpu=570;fs1=280;fsu=900
wp=[2*fp1/Fs,2*fpu/Fs];ws=[2*fs1/Fs,2*fsu/Fs];rp=0.1;rs=60;
[N,wp]=ellipord(wp,ws,rp,rs);
[B,A]=ellip(N,rp,rs,wp);
y2t=filter(B,A,st);
figure(3);subplot(2,1,1);
[H,W]=freqz(B,A,1000);
m=abs(H);
plot(W/pi,20*log10(m));grid on;
xlabel('\\omega/\\pi');ylabel('-A(dB)')
axis([0,1,-110,5]);title('损耗函数曲线');
figure(3);subplot(2,1,2);
n=0:length(y2t)-1;t=n*T;
plot(t,y2t);grid on;
xlabel('t/s');ylabel('yn');
axis([0,t(end),min(y2t),1.2*max(y2t)]);title('输出波形');
clear all;close all;
Fs=10000; T=1/Fs;
st=mstg;
fp=880;fs=550;
wp=2*fp/Fs;ws=2*fs/Fs;rp=0.1;rs=60;
[N,wp]=ellipord(wp,ws,rp,rs);
[B,A]=ellip(N,rp,rs,wp,'high');
y3t=filter(B,A,st);
figure(4);subplot(2,1,1);
[H,W]=freqz(B,A,1000);
m=abs(H);
plot(W/pi,20*log10(m));grid on;
xlabel('\\omega/\\pi');ylabel('-A(dB)')
axis([0,1,-100,5]);title('损耗函数曲线');
figure(4);subplot(2,1,2);
n=0:length(y3t)-1;t=n*T;
plot(t,y3t);grid on;
xlabel('t/s');ylabel('yn');
axis([0,t(end),min(y3t),1.2*max(y3t)]);title('输出波形');
信号产生函数xtg清单
function xt=xtg(N)
N=2000;Fs=1000;T=1/Fs;Tp=N*T;
t=0:T:(N-1)*T;
fc=Fs/10;f0=fc/10;
mt=cos(2*pi*f0*t);
ct=cos(2*pi*fc*t);
xt=mt.*ct;
nt=2*rand(1,N)-1;
fp=150;fs=200;Rp=0.1;As=70;
fb=[fp,fs];m=[0,1];
dev=[10^(-As/20),(10^(Rp/20)-1)/(10^(Rp/20)+1)];
[n,fo,mo,W]=remezord(fb,m,dev,Fs);
hn=remez(n,fo,mo,W);
yt=filter(hn,1,10*nt);
xt=xt+yt;
fst=fft(xt,N);k=0:N-1;f=k/Tp;
subplot(3,1,1);plot(t,xt);grid;xlabel('t/s');ylabel('x(t)');
axis([0,Tp/5,min(xt),max(xt)]);title('(a)信号加噪声波形')
subplot(3,1,3);plot(f,abs(fst)/max(abs(fst)));grid;title('(b)信号加噪声的频谱')
axis([0,Fs/2,0,1.2]);xlabel('f/Hz');ylabel('幅度')
滤波器的设计
function y = myplot(B,A);
[HK,WK] = freqz(B,A,1024);
plot(WK/pi,20*log10(abs(HK)));grid on;
xlabel('wk/pi');ylabel('幅度/db');
axis([0,1.0,-100,5]);
function tplot(xn,T,yn)
n = 0:length(xn)-1;
t=n*T;
plot(t,xn);
xlabel('t/s');ylabel(yn);
axis([0,t(end),min(xn),max(xn)]);
clear all;close all;
%指标给定
N = 1000;xt = xtg(N);
fp = 120;fs = 150;Rp = 0.2; As = 60;Fs = 1000;T = 1 / Fs;
%窗函数法设计低通滤波器
wc = (fp + fs) / Fs;
B = 2 * pi * (fs - fp) / Fs;
Nb = ceil(11 * pi / B);%blackman
hn = fir1(Nb - 1,wc, blackman(Nb));
Hw = abs(fft(hn, 1024));
ywt = fftfilt(hn, xt, N);
subplot(2,1,1);
myplot(hn,1);
%title('(a) 窗函数法低通滤波器幅频特性');
title('(a) 低通滤波器幅频特性');
subplot(2,1,2);
yt = 'y_w(t)';
tplot(ywt,T,yt);
title('(b) 滤除噪声后的波形');
%等纹波最佳逼近法设计低通滤波器
fb = [fp, fs];m = [1, 0];
dev = [(10^(Rp/20) - 1)/(10^(Rp/20) + 1), 10^(-As/20)];
[Ne, fo, mo, W] = remezord(fb, m, dev, Fs);
hn = remez(Ne, fo, mo, W);
Hw = abs(fft(hn, 1024));
yet = fftfilt(hn, xt, N);
%绘图
%freqz(hn,1,1024,Fs);
%title('Lowpass Filter Designed to Specifications');
figure(2);
subplot(2,1,1);
myplot(hn,1);
%title('(a) 等纹波最佳逼近法');
title('(a) 低通滤波器幅频特性');
subplot(2,1,2);
yt = 'y_e(t)';
tplot(ywt,T,yt);
title('(b) 滤除噪声后的波形');
修改信号后的st的时域波形和幅频特性曲线下载本文
