探究插值抗镜像 FIR 低通滤波器的参数对插值结果的影响

---

前言

实验结合了插值抗镜像 FIR 低通滤波器的整个设计过程，通过控制滤波器参数来具体分析其参数对插值后结果的影响。

---

一、实验目的

(1) 设计线性相位的 FIR 滤波器实现插值抗镜像滤波器。
(2) 分析 FIR 阶数和截止频率等参数对插值信号的影响。

二、实验过程

1. 原理分析与程序设计（框图与描述）

插值是提高采样率以增加数据的过程，在实际实现过程中采取先插零再滤波的方法来达到这个目的。流程如下，信号 x(n)先通过插零，即在x(n)中每个点后面插入（I-1）个零数据点，插零会导致镜像频率的出现，即频谱会以采样频率的整数倍为中心左右延拓（N ∙ fs ± f），


故要通过滤波器滤除这些多余频率，插零后信号再通过设计的 FIR 滤波器滤波后得到 y(n)，这时采样率变为原来的 I 倍。

对于 FIR 滤波器，本实验采用了等波纹最佳逼近法来设计实验所需的 FIR滤波器，在 matlab 中采用 remezod 和 remez 函数来快速设计采用等波纹最佳逼近法的 FIR 滤波器。
其使用方法如下：

fb = [fp,fst]; m= [1,0]; % 确定 remezord 函数所需参数 f,m
rip = [(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)]; % 参照公式计算通带和阻
带波纹幅度
% 确定 remez 函数所需参数，其中 N 为滤波器阶数, Fs 为滤波器采样频率，
可更改
[N,fo,mo,W] = remezord(fb,m,rip,Fs); % 估算滤波器参数
hn = remez(N,fo,mo,W); % 调用 remez 函数构成滤波器
Hw = abs(fft(hn,1024)); % 求设计的滤波器频率特性

其中，f 提供通带和阻带截止频率，m 提供与 f 对应幅度特性，f 和 m 一起给 remezord 函数提供希望逼近的幅度特性曲线；N 为滤波器阶数，remezord 函数返回的 N 是所设定参数下估计的最小滤波器阶数，在探究 FIR 阶数对插值信号的影响时，可以更改 remez 输入参数中的 N。

2.程序设计（框图与描述）

代码流程为先生成一个简单的正弦采样信号（方便滤波器参数选择），然后对其进行插零，插零后进行 FIR 抗镜像滤波，滤波完成后绘制插值后的信号波形图和频谱。

主代码为：

clear;close all;clc;
% 设置参数
I = 5; % 插值因子

% 选择1：生成信号
Tp = 0.5; % 持续时间
fs = 2000; % 采样频率
Fs = fs*I; % 插值后频率
f = 100; % 信号频率
t = 0:1/fs:Tp-1/fs; % 时间
y = sin(2*pi*f*t);
% 绘制原信号幅度谱
figure;
subplot(2,1,1);
stem(t,y);
axis([0 0.01 -1 1]);
xlabel('时间 (s)');
ylabel('幅度');

% 绘制频谱图
subplot(2,1,2);
freqStem(y,fs);

% 插零
y_len = length(y);
y_1 = [];
for ii = [1:y_len]
    y_1 = [y_1,y(ii),zeros(1,I-1)];
end
% 绘制插零后幅度谱
figure;
subplot(2,1,1);
t_1 = 0:1/Fs:Tp-1/Fs;
stem(t_1,y_1);
axis([0 0.01 -1 1]);
xlabel('时间 (s)');
ylabel('幅度');
% 绘制频谱图
subplot(2,1,2);
freqStem(y_1,Fs); % 信号频率会以fs为的整数倍为中心左右延拓形成镜像频率

% 设定指标
fp = 1000; fst = 1500;  % 通带和阻带截止频率(可更改)
Rp = 0.2; As = 60;  % 通带最大衰减，阻带最小衰减

% 1.用窗函数法设计FIR滤波器
% wc = (fp+fst)/Fs;     % 理想低通滤波器截止频率(关于pi归一化）
% B = 2*pi*(fst-fp)/Fs;  % 过渡带宽度指标
% % Nb = ceil(11*pi/B);    % blackman窗的长度N
% Nb = 128;
% hn = fir1(Nb-1,wc,blackman(Nb));
% Hw = abs(fft(hn,1024));	% 求设计的滤波器频率特性

% 2.用等波纹最佳逼近法设计FIR滤波器
fb = [fp,fst]; m= [1,0]; 		% 确定remezord函数所需参数f,m
rip = [(10^(Rp/20)-1)/(10^(Rp/20)+1),10^(-As/20)]; % 参照公式计算通带和阻带波纹幅度
% 确定remez函数所需参数，其中N为滤波器阶数,Fs为滤波器采样频率，可更改
[N,fo,mo,W] = remezord(fb,m,rip,Fs);	
hn = remez(N,fo,mo,W);		% 调用remez函数构成滤波器
Hw = abs(fft(hn,1024));		% 求设计的滤波器频率特性

% 绘制滤波器频率响应曲线
f = [0:1023]*Fs/1024;
figure;
subplot(3,1,1);
plot(f,20*log10(Hw/max(Hw)));grid;
axis([0,Fs/2,-80,10]);
xlabel('f/Hz');ylabel('幅度');
str=sprintf('fst=%d',fst);
title('等波纹法低通滤波器幅频特性（N=50，fst=1.5kHz）');

% 调用函数fftfilt对y_1滤波
filteredSignal = I*fftfilt(hn,y_1);       

% 绘制滤波后幅度谱
subplot(3,1,2);
stem(t_1,filteredSignal);
axis([0 0.015 1.1*min(filteredSignal) 1.1*max(filteredSignal)]);
xlabel('时间 (s)');
ylabel('幅度');
% 绘制频谱图
subplot(3,1,3);
freqStem(filteredSignal,fs*I);

在本实验中，为了绘制频谱图方便，我定义了一个函数： freqStem(x,fs,n)

function freqStem(x,fs,n)
if nargin > 2
    str=sprintf('n=%d',n);
    title(str + '的频谱 ');
end
N = length(x);
X = fftshift(fft(x)); % 进行FFT操作并移位
freqAxis = (-fs/2) + ((0:length(X)-1) * (fs/length(X)));
plot(freqAxis, 2*abs(X)/N);
axis([-fs/2 fs/2 0 1.2*max(2*abs(X)/N)]);
xlabel('频率 (Hz)');
ylabel('幅度');

可以直接绘制信号频谱。其中，x 指输入信号，fs 指信号采样率，n 指信号名（可空）。

3.实验结果与分析

在该实验中我生成了一个采样频率为 2kHz，信号频率为 100Hz 的采样点数
为 1k 的正弦采样信号信号，其波形图（时间显示范围进行了设置）和频谱如下。

本实验取 I =5，信号插零后（即每个点后插入 4 个零）波形图和频谱如下，可见其频谱以采样频率的整数倍为中心左右延拓。经插零后，采样点数变为了5k。同时也可看出，频率幅度变为了原来的 1/I，即 1/5。

为了滤除这些多余频率，FIR 滤波器采用了通带截止频率为 1kHz，阻带截止频率为 1.5kHz，N 的估算最小值为 50。经滤波后信号波形图（这里方便观察，对幅度进行了归一化，正常情况滤波后信号幅度会变为原来的 1/I，即 1/5）和频谱如下，除红框标记（由于运行速率原因，一开始输入的信号个数小于滤波器的所需要的数据，为无效区）外，波形与原始信号一致。在 matlab工作区可见其长度即采样点数为 5k 点，信号频率保持不变，采样率变为了原来的 I 倍，即 10kHz。

(1) 阶数 N 对结果影响：
当滤波器 fst 保持不变，阶数 N 增大为 100，无效区也会变大，约为原来 2倍。波形依旧保持不变。

当滤波器 fst 保持不变，阶数 N 减小为 10 时，无效区虽然减小，但滤波器滤波不完全，出现了很多小的频率分量，使波形不完整。

(2) 阻带截止频率 fst 对结果影响
当滤波器阶数 N 保持为 50 不变，fst 增大到 1.9kHz 以上（这里用 2.5kHz）时，其通带中会含有需滤除的分量，使滤波后信号频谱会出现其它频率分量，没有到达滤除镜像频率的目的，使结果波形不符合要求。

---

总结

通过这部分实验，可以了解通过插值来提高信号采样率的基本概念，理解如何设计相对应的FIR 滤波器滤除这些镜像频率，在探究 FIR 滤波器参数对实验结果的影响过程中可以更好地理解 FIR 滤波器对参数的选择。