lowpassfir.m

交互式低通滤波器

function y=lowpassfir(sn,Fc)
%LOWPASS Returns a discrete-time filter object.
%
% M-File generated by MATLAB(R) 7.1 and the Signal Processing Toolbox 6.4.
%
% Generated on: 30-May-2008 16:33:05
%

% Equiripple Lowpass filter designed using the FIRPM function.

% All frequency values are in Hz.
Fs = 5*Fc;               % Sampling Frequency
Fpass = Fc-20;           % Passband Frequency
Fstop = Fc+20;           % Stopband Frequency
Dpass = 0.057501127785;  % Passband Ripple
Dstop = 0.0001;          % Stopband Attenuation
dens  = 20;              % Density Factor

% Calculate the order from the parameters using FIRPMORD.
[N, Fo, Ao, W] = firpmord([Fpass, Fstop]/(Fs/2), [1 0], [Dpass, Dstop]);
% Calculate the coefficients using the FIRPM function.
b  = firpm(N, Fo, Ao, W, {dens});
Hd = dfilt.dffir(b);
% [EOF]
N=length(sn);                    %检测带噪声信号的长度
t=0:N-1;
y=filter2(Hd.Numerator,sn);      %对信号sn滤波，y是滤波后的信号
figure
plot(t,sn,'r',t,y,'k',t,y-sn,'b');
legend('带噪声信号','滤波后信号','噪声');
%****对比去噪声前后信号的频谱图*****************
M=pow2(nextpow2(N));              %求最小的大于N的2的幂数
f = Fs*(0:M/2-1)/M;              %把序列转化成频率
nf=length(f);  
%带噪声信号
Msn=fft(sn,M);                     %FFT变换得到复数形式的序列；
Pyys = Msn.* conj(Msn)/M;         %序列中每个元素和它的共轭相乘得到相应的幅值
ps = unwrap(angle(Msn));          %求序列每个复数元素的主值角
%去噪声后信号
My=fft(y,M);
Pyyy = My.* conj(My)/M;
py = unwrap(angle(My));
%画图幅度特性
figure
subplot(1,2,1);
plot(f,Pyys(1:nf));
xlabel('frequency (Hz)')
legend('带噪声信号');
subplot(1,2,2);
plot(f,Pyyy(1:nf));
xlabel('frequency (Hz)')
legend('滤波后信号');
title('幅度特性');
%画图相位特性
figure
subplot(1,2,1);
plot(f,ps(1:nf));
xlabel('frequency (Hz)')
legend('带噪声信号');
subplot(1,2,2);
plot(f,py(1:nf));
xlabel('frequency (Hz)')
legend('滤波后信号');
title('相位特性');