matlab fir低通滤波器 firpm_信号叠加、分析及信号滤波频谱—FIR,IIR

一、实验目的

(1)综合应用信号频谱分析和数字滤波器设计的知识，实现信号的滤波。

(2)比较分析干扰信号加入前、后的频谱有何区别，分析噪声的影响。

(3)加深理解信号时域和频域分析的物理概念，认识不同滤波器的特性和适用范围。

(4)掌握设计IIR、FIR 数字滤波器的原理和特性，通过比较对两种滤波器的滤波作用有直观认识。

二、实验原理

(1) FIR数字滤波器的设计原理：

FIR数字滤波器总是稳定的系统，且可以设计成具有线性相位。M阶FIR数字滤波器的系统函数为

系统的单位脉冲响应h[k]是长度为M+1的有限长因果序列。当满足

h[k] =±h[M -k]的对称条件时，该FIR数字滤波器具有线性相位。FIR数字滤波器设计方法主要有窗口法、频率取样法、优化设计法。

MATLAB 中常用FIR数字滤波器设计函数：

fir1 窗函数法FIR数字滤波器设计：低通、高通、带通、带阻、多频带滤波器。

fir2 频率抽样法FIR数字滤波器设计：任意频率响应。

firls FIR数字滤波器设计：指定频率响应。

(2) IIR 数字滤波器的设计原理：

IIR数字滤波器一般为线性时不变的因果离散系统，N阶IIR数字滤波器的系统函数可以表达为z-1的有理多项式，即

分母系数中{

；i=1,2,…,N}至少有一个非零。对于因果IIR数字滤波器，满足M≤ N。

IIR数字滤波器的设计主要通过成熟的模拟滤波器设计方法来实现。首先在频域将数字滤波器设计指标转化为模拟低通滤波器设计指标，将任意模拟滤波器转换为原型模拟低通滤波器设计指标，根据模拟低通滤波器设计指标设计出模拟低通滤波器Hlp(s)，由Hlp(s)经过相应的复频域变换得到H(s)，由H(s)经过脉冲响应不变法或双线性变换法得到所需的IIR数字滤波器H(z)，由此可见，IIR数字滤波器的设计的重要环节是模拟低通滤波器的设计，涉及低通模拟滤波器的方法有巴特沃斯，切比雪夫和椭圆等滤波器设计方法。

用 MATLAB 设计IIR滤波器的常用函数：

IIR滤波器阶数选择：

Buttord-------巴特沃斯（Butterworth）滤波器阶数选择。

Cheblord-----切比雪夫（Chebyshev）I 型滤波器阶数选择。

Cheb2ord----切比雪夫（Chebyshev）II 型滤波器阶数选择。

Ellpord-------椭圆（Elliptic）滤波器阶数选择。

IIR滤波器设计：

Butter-------巴特沃斯（Butterworth）滤波器设计。

cheby1------切比雪夫（Chebyshev）I 型滤波器设计。

cheby2------切比雪夫（Chebyshev）II 型滤波器设计。

maxflat------通用的巴特沃斯（Butterworth）低通滤波器设计。

三、实验步骤

1利用 wavread 函数读入一段音乐信号(.wav)，利用 FFT 分析该信号的频谱。

[y,FS]=audioread('C:UsersDesktopwind.mp3');
music=y(:,1);   %由于原MP3音频为双音轨，除去一条
N=length(music);    %获取原音频的采样点数
w=linspace(0,2*pi,N);   %设置数字角频率范围
MUSIC=fft(music);   %求原音乐信号频谱
subplot(2,2,1)
pMUSIC=abs(MUSIC);  %绘制原信号幅度谱
plot(w,pMUSIC);
title('原始音乐频谱');
xlabel('数字角频率/rad');
ylabel('幅度模值');
audiowrite('wind1.wav',music,FS);  %保存除去一条音轨的文件
T=1/FS; %开始描述噪声信号
w1=2*pi*500*T;  %将模拟角频率转化为数字角频率
w2=2*pi*1000*T;
for n=1:N
    noise(n)=(sin(w1*n)+sin(w2*n))*0.1; %调整合适噪声大小，不要把原信号全部淹没
end
NoiseMusic=music+noise';    %将噪声加入原信号
subplot(2,2,2)
pnoise=abs(fft(noise)); %绘制噪声频谱
plot(w,pnoise);    %加abs画幅度谱
title('噪声频谱');
xlabel('数字角频率/rad');
ylabel('幅度模值');
subplot(2,2,3)
pNoiseMusic=abs(fft(NoiseMusic));
plot(w,pNoiseMusic);
title('加入噪声后音乐频谱');
xlabel('数字角频率/rad');
ylabel('幅度模值');
audiowrite('windnoise.wav',NoiseMusic,FS);  %保存加入噪声的音频
Rp=2;   %设置通带波纹为2dB
Rs=60;  %设置阻带衰减为60dB
%除去500Hz干扰信号
wp=[20,700]*2/FS;   %设置通带区间
ws=[499,501]*2/FS;  %设置阻带区间
[n,wn]=buttord(wp,ws,Rp,Rs);    %巴特沃斯滤波器阶数选择
[b,a]=butter(n,wn,'stop');  %设计带阻巴特沃斯滤波器
music1=filter(b,a,NoiseMusic); %求经过除去500Hz带阻后的信号
%除去1000Hz干扰信号
wp1=[900,1300]*2/FS;
ws1=[999,1001]*2/FS;
[n,wn]=buttord(wp1,ws1,Rp,Rs);
[b1,a1]=butter(n,wn,'stop');
music2=filter(b1,a1,music1);
%绘制除去500Hz和1000Hz干扰信号后的音乐频谱
subplot(2,2,4)
pmusic2=abs(fft(music2)); 
plot(w,pmusic2);
title('除去噪声后的频谱');
xlabel('数字角频率/rad');
ylabel('幅度模值');
audiowrite('windpure.wav',NoiseMusic,FS);   %保存除噪声后的信号

2在该段音乐信号中加入两个频率分别为 500Hz 和 1000Hz 的正弦型干扰信号，利用 FFT 观察该信号的频谱。

3根据其频谱，设计不同类型的滤波器（IIR，FIR）滤除干扰信号，从时域（通过播放）和频域（通过频谱）比较结果。

通过播放之后发现滤波效果是非常不错的，噪声已经基本被滤去。

上为两个文件，一个为滤波之前加入噪声的音乐文件windnoise.wav,后者为滤波之后的windpure.wav.

四、思考题

1.利用频率选择滤波器（FIR 滤波器、IIR 滤波器）进行信号去噪的基本思想及主要步骤是什么？

基本思想：利用频率选择滤波器的传输选择特性，对输入信号进行加工和变换，改变输入序列的频谱或信号波形，使需要的频率分量通过，抑制无用的信号分量。

主要步骤：首先根据输入信号的频谱特性及设计目标要求确定需要使用的滤波器类型，包括低通、高通、带通、带阻。其次，由需要的相位特性及幅度特性确定滤波器的设计方法，包括FIR（巴特沃斯、切比雪夫）、IIR(巴特沃斯、切比雪夫)，然后根据目标设计要求确定滤波器参数，需要的参数根据滤波器的设计方法确定，主要包括滤波器的阶次、频率要求、幅度响应要求，之后再利用设计好的滤波器对实际信号进行处理，并检验是否满足设计要求，并调整和优化。

2.如何根据含有噪声信号的频谱特性选择滤波器的类型和设计指标？

根据采集到的信号获得频谱图，由时域频域对应关系，确定需要滤除噪声的特性，然后从以下角度确定所需要的滤波器：首先是频响特性，IIR滤波器设计时不考虑相位特性，且通常相位都是非线性的，而FIR滤波器在满足幅频特性要求的同时，还能获得比较严格的线性相位特性，利用窗函数或者其他算法可以逼近更加任意的频响特性，因此性能更优越，适用范围更广。其次是稳定性问题，IIR滤波器设计时，极点必须在单位圆之内，而FIR滤波器极点在单位圆内，因此始终稳定。然后考虑滤波器的结构，IIR滤波器一般采用递归结构，存在有输出对输入的反，而IIR滤波器的阶次相对较低，运算次数少，存储单元少，FIR滤波器正好相反。最后考虑设计工作量，FIR无表可查，需要用到迭代法，计算量较大，而IIR滤波器相对简单，有现成的计算公式和数据表格可以查阅。

实验二 图像滤波实验

一、实验目的

设计一个9696像素的离散空间图像。假定图像为一个棋盘，由8*8的交替变化的黑白方块组成，如图。

（a） 用冲激响应为

的滤波器对图像进行逐行滤波，然后逐列滤波。用matlab的imagesc命令在屏幕上显示图像。

（b） 用冲激响应为

的滤波器对图像进行滤波，然后进行逐列滤波。用matlab的imagesc命令在屏幕上显示图像。

二、实验结果

N=96;   %设置分辨率
p=ones(N,N);
for m=1:96  %将题目中所述图像二值化
    if mod(floor((m-1)/12),2)==0
    for n=1:96
        if mod(floor((n-1)/12),2)==0
            p(m,n)=0;
        end
    end
    else
        for n=1:96
           if mod(floor((n-1)/12),2)==1
            p(m,n)=0;
           end
        end
    end
end
colormap gray
subplot(1,3,1)
imshow(p);  %绘制原图像
title('原图像');
for n=1:96
    h1(n)=0.2*0.8^n;    %求(a)中的冲激响应
end
h2(1)=1-0.2*0.8;    %求(b)中的冲激响应
for n=2:96
   h2(n)=-0.2*0.8^n;
end
%接下来求变化后图像
for i=1:96
   p1(i,:)=ifft(fft(p(i,:)).*fft(h1)); %求经过h1逐行滤波后的图像
   p2(i,:)=ifft(fft(p(i,:)).*fft(h2)); %求经过h2逐行滤波后的图像
end
h11=h1';
h21=h2';
for i=1:96
    p1(:,i)=ifft(fft(p1(:,i)).*fft(h11)); %求经过h1逐列滤波后的图像
    p2(:,i)=ifft(fft(p2(:,i)).*fft(h21)); %求经过h2逐列滤波后的图像
end
subplot(1,3,2)
imshow(p1);
title('经过(a)所示滤波器图样');
subplot(1,3,3)
imshow(p2);
title('经过(b)所示滤波器图样');