2019-12-19 16:45:03 hopeless123 阅读数 221

最近需要用到Fir滤波器,在网上也看了不少资料,发现一个稍微能用的(https://blog.csdn.net/BIGFatming/article/details/92386914),主要代码也是直接copy的,但是在使用过程中发现,三角窗的实现好像不对,而且只实现了低通的,根据该文章内容,我自己重写了一个,包括三角窗函数的实现,添加了高通、带通、带阻的实现。环境Visual Studio2015,c#。

1、如下图(最下面的曲线),使用网友的三角函数实现不对;

2、据此更改了相关代码;

public class Bartlett
    {
        public int N = 0;
        public Bartlett(double Wp, double Ws)
        {
            int i;
            double n = (2.1 * 2 * Math.PI) / (Ws - Wp);
            for (i = 0; i < n; i++) ;
            N = i;
        }

        public double[] GetWin()
        {
            int n;
            double[] wd = new double[N];
            for (n = 0; n < N; n++)
            {
                if (n <= ((N - 1) / 2))
                {
                    wd[n] = (double)(2 * (double)n / ((double)N - 1));
                }
                if (n > ((N - 1) / 2))
                {
                    wd[n] = 2 * ((double)N - 1 - n) / ((double)N - 1);
                }
            }
            return wd;
        }
}

3、下图为Blackman和三角函数的h(n);

4、下图为Hanning和Hamming窗的h(n);

5、由于看到的资料只提供了低通的实现方式,根据《数字信号处理理论、算法与实现(第三版)》实现了对应的高通、带通、带阻的实现;

6、高通

public double[] GetHigh()
        {
            double[] hd = new double[N];
            for (int n = 0; n < N; n++)
            {
                double numerator = Math.Sin(Math.PI * (n - alpha)) - Math.Sin((n - alpha) * Wc);
                double denominator = Math.PI * (n - alpha);
                if (n == alpha)
                {
                    hd[n] = (Math.PI - Wc) / Math.PI;
                }
                else
                {
                    hd[n] = numerator / denominator;
                }
            }
            return hd;
        }

7、带通

public double[] GetDaiTong(double Wl, double Wh)
        {
            double[] hd = new double[N];
            for (int n = 0; n < N; n++)
            {
                double numerator = Math.Sin((n - alpha) * Wh) - Math.Sin((n - alpha) * Wl);
                double denominator = Math.PI * (n - alpha);
                if (n == alpha)
                {
                    hd[n] = (Wh - Wl) / Math.PI;
                }
                else
                {
                    hd[n] = numerator / denominator;
                }
            }
            return hd;
        }

8、带阻,这里设计时,对于求n=alpha时的hd[n]值不是很确定,上图中用红笔注释的地方(为下面代码实现的内容)。

注:能力有限,在分析原始算法时,我是无法求解n=alpha的值或表达式(有人知道原始算法如何求解的话请留言,3Q),故带阻的实现可能在n=alpha点时是错的。【后期再进行研究,及时更新】。

public double[] GetDaiZu(double Wl, double Wh)
        {
            double[] hd = new double[N];
            for (int n = 0; n < N; n++)
            {
                double numerator = Math.Sin((n - alpha) * Wl) + Math.Sin(Math.PI * (n - alpha)) - Math.Sin((n - alpha) * Wh);
                double denominator = Math.PI * ((double)n - alpha);
                if (n == alpha)
                {
                    hd[n] = (Math.PI - Wh) / Math.PI + Wl / Math.PI;
                }
                else
                {
                    hd[n] = numerator / denominator;
                }
            }
            return hd;
        }

9、测试:(测试代码均使用的hanning窗,测试低通、高通、带通与带阻功能)

9.1、原始信号参数为:Fs=1000,采样点数:2048,原始数据包含两个频率点数据5Hz和20Hz。Matlab代码:

f1=5;f2=20;f3=15;
t=0:0.001:2.049;
y=2*cos(2*pi*f1*t)+2*sin(2*pi*f2*t);
plot(t,y);
fid=fopen('test1.txt','wt');
fprintf(fid,'%d,',y);
fclose(fid);

9.2、原始信号及频谱,

9.3、低通滤波,Wp=0.02 * Math.PI,Ws=0.03 * Math.PI(分别对应Wp=10Hz,Ws=15Hz)

注:由于给入的频率为归一化后的,需要计算((Fs/2)*Wp=fp,截止频率为fc【Wc为归一化频率】);相位延迟:(1/Fs)*(N-1)/2;滤波器阶数:N-1(N为窗长度)【这里不确定,网上查到有这个说法,有大神知道的话可以探讨一下】。

9.4、高通滤波,参数与低通一致,

9.5、带通滤波,带通频率设置为【0.03 * Math.PI, 0.05 * Math.PI】(即为15Hz~25Hz,通带)

9.6、带阻滤波,带通频率设置为【0.03 * Math.PI, 0.05 * Math.PI】(即为15Hz~25Hz,阻带)

10、总结:

1)基本能实现滤波功能,但是有很多参数以及算法需要再研究优化,仅提供一种思路;(未涉及到FIR的4种形式)

2)从效果来看,高通与带阻滤波后会出现很多高频的干扰(频段),而低通与带通没有,图中能看出,还需要进一步研究;

3)本文采用的是时域卷积的处理方式,其实也可以用差分方程来做,自己赖还没弄,处理完后都有数据点增多以及相移(线性)的问题,需要想办法优化;

4)整个测试代码https://download.csdn.net/download/hopeless123/12040952。

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

        

 

 

 

 

2019-11-16 01:13:26 qq_17119267 阅读数 666

[Matlab]FIR滤波器设计:(基本窗函数FIR滤波器设计)

​ IIR滤波器主要设计方法先设计一个模拟低通滤波器,然后把它转化为形式上的数字滤波器。但对于FIR滤波器来说,设计方法的关键要求之一就是保证线性相位条件。而IIR滤波器的设计方法中只对幅值特性进行设计,因此无法保证相位。所以FIR滤波器的设计需要采用完全不同的方法。FIR滤波器的设计方法主要有窗函数法、频率采样法、切比雪夫逼近法等。

​ 由于理想滤波器在边界频率处不连续,故时域信号hd(n)h_d(n)一定是无限时宽的,也是非因果的。所以理想低通滤波器是无法实现的。如果实现一个具有理想线性相位特性的滤波器,其幅值特性只能采用逼近理想幅频特性的方法实现。如果对时域信号hd(n)h_d(n)进行截取,并保证截取过程中序列保持对称,而且截取长度为N,则对称点为α=12(N1)α=\frac{1}{2}*(N-1)。截取后序列为h(n)h(n),侧h(n)h(n)可用下式子表示:
h(n)=hd(n)w(n) h(n) = h_d(n)*w(n)
式中,w(n)w(n)为截取函数,又称为穿函数。如果窗函数为矩形序列,则称之为矩形窗。窗函数有多种形式,为保证加窗后系统的线性相位特性,必须保证加窗后序列关于α=12(N1)α=\frac{1}{2}*(N-1)点对称。常用的窗函数有矩形窗、汉宁窗、汉明窗、布莱克曼窗、凯塞窗。窗函数设计法的基本思想是用一个长度为N的序列h(n)h(n)代替hd(n)h_d(n)作为实际设计的滤波器的单位脉冲响应。这种设计法成为窗函数设计法。显然在保证h(n)h(n)对称性的前提下,窗函数长度N越长,则h(n)h(n)越接近hd(n)h_d(n)。但是误差是肯定存在的,这种误差成为截断误差。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-A0flZjTw-1573836815510)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\20170402112228645.jpg)]

​ 要确定如何设计一个FIR滤波器,首先得对加窗后的理想滤波器特性变化进行分析,并研究减少由截断引起的误差的途径,从而提出更好的滤波器设计方案。对于调整窗口长度可以有效地控制过渡带的宽度,但减少带内波动以及加大阻带衰减没有作用。所以必须挑选最为合适的窗函数对理想滤波器进行截取。下面简单介绍几种窗函数。一个实际的滤波器的单位脉冲响应可表示为:
h(n)=hd(n)w(n) h(n) = h_d(n)*w(n)
几种常见窗函数(只给出了窗函数的定义和幅度特性):
W(ejw)=W(w)(ejαw) W(e^{j*w})=W(w)(e^{-jαw})

矩形窗FIR滤波器设计:

矩形窗的窗函数为:
wR(n)=RN(n) w_R(n)=R_N(n)
幅度函数为:
RN(w)=sin(wN/2)sin(w/2) R_N(w) = \frac{sin(wN/2)}{sin(w/2)}
它的主瓣宽度为4π/N4\pi/N,第一瓣比主瓣地13dB.

在Matlab中,实现矩形窗的函数为boxcar和recttwin ,其调用格式如下:

w=boxcar(N);
w=recttwin(N);
%%%显示窗函数的GUI工具
n = 60;
w = rectwin(n);
wvtool(w)%显示窗函数的GUI工具
%还提供了显示窗函数的GUI工具,如wvtool可以显示用来显示窗的形状和频域图形,wintool可以打开窗设计和分析工具,如运行
wvtool(hamming(64),hann(64),gausswin(64))
%%可以对比汉明窗、汉宁窗和高斯窗

其中,N是窗函数的长度,返回值w是一个N阶的向量,它的元素有窗函数的值组成。其中w=boxcar(N)等价于w=ones(N,1)。

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-eknDhIxt-1573836815512)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\20170402112307695.png)]

案例分析:

利用矩形窗设计FIR带阻滤波器,主要参数如下:

给定抽样频率为Ωs=2pi1.5104(rad/sec)Ωs=2*pi*1.5*10^4(rad/sec),

通带截至频率为Ωp1=2pi0.75103(rad/sec),Ωp2=2pi6103(rad/sec)Ωp1=2*pi*0.75*10^3(rad/sec),Ωp2=2*pi*6*10^3(rad/sec)
阻带截至频率为Ωst1=2pi2.25103(rad/sec),Ωst2=2pi1.5103(rad/sec)Ωst1=2*pi*2.25*10^3(rad/sec),Ωst2=2*pi*1.5*10^3(rad/sec)
阻带衰减δ2>=50dB{\delta}_2 >=50dB

%%%%调用子程序1:
function hd=ideal_bs(Wcl,Wch,m); 
alpha=(m-1)/2; 
n=[0:1:(m-1)];
m=n-alpha+eps; 
hd=[sin(m*pi)+sin(Wcl*m)-sin(Wch*m)]./(pi*m)
%%%%调用子程序2:
function[db,mag,pha,w]=freqz_m2(b,a)
[H,w]=freqz(b,a,1000,'whole');
H=(H(1:1:501))'; w=(w(1:1:501))';
mag=abs(H);
db=20*log10((mag+eps)/max(mag));
pha=angle(H);
%%%%运行MATLAB源代码如下:
clear all;
Wph=3*pi*6.25/15;%通带频率
Wpl=3*pi/15;
Wsl=3*pi*2.5/15;%阻带频率
Wsh=3*pi*4.75/15;
tr_width=min((Wsl-Wpl),(Wph-Wsh));%%过渡带带宽
%过渡带宽度
N=ceil(4*pi/tr_width);					%滤波器长度
n=0:1:N-1;
Wcl=(Wsl+Wpl)/2;						%理想滤波器的截止频率
Wch=(Wsh+Wph)/2;
hd=ideal_bs(Wcl,Wch,N);				%理想滤波器的单位冲击响应
w_ham=(boxcar(N))';
string=['矩形窗','N=',num2str(N)];
h=hd.*w_ham;						%截取取得实际的单位脉冲响应
[db,mag,pha,w]=freqz_m2(h,[1]);
%计算实际滤波器的幅度响应
delta_w=2*pi/1000;
subplot(241);
stem(n,hd);
title('理想脉冲响应hd(n)')
axis([-1,N,-0.5,0.8]);
xlabel('n');ylabel('hd(n)');
grid on
subplot(242);
stem(n,w_ham);
axis([-1,N,0,1.1]);
xlabel('n');ylabel('w(n)');
text(1.5,1.3,string);
grid on
subplot(243);
stem(n,h);title('实际脉冲响应h(n)');
axis([0,N,-1.4,1.4]);
xlabel('n');ylabel('h(n)');
grid on
subplot(244);
plot(w,pha);title('相频特性');
axis([0,3.15,-4,4]);
xlabel('频率(rad)');ylabel('相位(Φ)');
grid on
subplot(245);
plot(w/pi,db);title('幅度特性(dB)');
axis([0,1,-80,10]);
xlabel('频率(pi)');ylabel('分贝数');
grid on
subplot(246);
plot(w,mag);title('频率特性')
axis([0,3,0,2]);
xlabel('频率(rad)');ylabel('幅值');
grid on
fs=15000;
t=(0:100)/fs;
x=cos(2*pi*t*750)+cos(2*pi*t*3000)+cos(2*pi*t*6100);
q=filter(h,1,x);
[a,f1]=freqz(x);
f1=f1/pi*fs/2;
[b,f2]=freqz(q);
f2=f2/pi*fs/2;
subplot(247);
plot(f1,abs(a));
title('输入波形频谱图');
xlabel('频率');ylabel('幅度')
grid on
subplot(248);
plot(f2,abs(b));
title('输出波形频谱图');
xlabel('频率');ylabel('幅度')
grid on

在这里插入图片描述

汉宁窗FIR滤波器设计:

汉宁窗(hanning window)又称升余弦窗,窗函数为:
wHn(n)=0.5[1cos(2πnN1)]RN(n) w_{Hn}(n) = 0.5[1-cos(\frac{2\pi*n}{N-1})]*R_N(n)
幅值函数为:
WHn(w)=0.5WR(w)+0.25[WR(w2πN1)+WR(w+2πN1)] W_{Hn}(w) = 0.5W_R(w) + 0.25[W_R(w-\frac{2\pi}{N-1})+W_R(w+\frac{2\pi}{N-1})]
汉宁窗幅度函数由3部分相加而成,其结果是使主瓣集中了更多能量,而旁瓣3部分相加时相互抵消而变小,其代价是主瓣宽度增加到8π/N8\pi/N。第一瓣比主瓣低31dB,阻带衰减加大。

在Matlab中,实现汉宁窗的函数为hanning和barthannwin ,其调用格式如下:

w=hanning(N)
w=barthannwin(N)

案例1:绘制50个点的汉宁窗。

N=49;n=1:N;
wdhn=hanning(N);	
figure(3);
stem(n,wdhn,'.');
grid on
axis([0,N,0,1.1]);
title('50点汉宁窗');
ylabel('W(n)');
xlabel('n');
title('50点汉宁窗');

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-KhfUUgrt-1573836815515)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanning_50point.bmp)]

案例2:已知连续信号为x(t)=cos(2πf1t)+0.15cos(2πf2t)x(t)=cos(2{\pi}f_1t) +0.15cos(2{\pi}f_2t),其中f1=100Hz,f2=150Hzf_1=100Hz,f_2=150Hz。若抽样频率fsam=600Hzf_sam=600Hzd对信号进行抽样,利用不同宽度N的矩形截断该序列,N取40,观察不同的窗对普分析结果的影响。

N=40; 
L=512;
f1=100;f2=150;fs=600;
ws=2*pi*fs;
t=(0:N-1)*(1/fs);
x=cos(2*pi*f1*t)+0.25*sin (2*pi*f2*t);
wh=boxcar(N)';
x=x.*wh;
subplot(221);stem(t,x);
title('加矩形窗时域图');
xlabel('n');ylabel('h(n)')
grid on
W=fft(x,L);
f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi); 
subplot(222);
plot(f,abs(fftshift(W)))
title('加矩形窗频域图');
xlabel('频率');ylabel('幅度')
grid on
x=cos(2*pi*f1*t)+0.15*cos(2*pi*f2*t); 
wh=hanning(N)';
x=x.*wh;
subplot(223);stem(t,x);
title('加汉宁窗时域图');
xlabel('n');ylabel('h(n)')
grid on
W=fft(x,L);
f=((-L/2:L/2-1)*(2*pi/L)*fs)/(2*pi);
subplot(224);
plot(f,abs(fftshift(W)))
title('加汉宁窗频域图');
xlabel('频率');ylabel('幅度')
grid on

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-wc7oZ6mR-1573836815515)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\rectange_6.8.bmp)]

用汉宁窗对谐波信号进行分析:

clear; 
% 原始数据:直流:0V; 基波:49.5Hz,100V,10deg; HR2:0.5V,40deg;
hr0=0;f1=50.1; 
hr(1)=25*sqrt(2);deg(1)=10; 
hr(2)=0;deg(2)=0; 
hr(3)=1.755*sqrt(2);deg(3)=40; 
hr(4)=0;deg(4)=0; 
hr(5)=0.885*sqrt(2);deg(5)=70; 
hr(6)=0;deg(6)=0; 
hr(7)=1.125;deg(7)=110; 
M=7;f=[1:M]*f1;							%设定频率 
% 采样 
fs=10000; 
N=2048;									% 约10个周期 
T=1/fs; 
n=[0:N-1];t=n*T; 
x=zeros(size(t)); 
for k=1:M 
    x=x+hr(k)*cos(2*pi*f(k)*t+deg(k)*pi/180); 
end 
%分析: 
w=0.5-0.5*cos(2*pi*n/N);
Xk=fft(x.*w); 
amp=abs(Xk(1:N/2))/N*2;						%幅频 
pha=angle(Xk(1:N/2))/pi*180;					%相频 
for k=1:N/2 
    if(amp(k)<0.01) pha(k)=0;  %当谐波<10mV时,其相位=0 
    end 
    if(pha(k)<0) pha(k)=pha(k)+360;%调整到0-360度 
    end 
end 
fmin=fs/N; 
xaxis=fmin*n(1:N/2);
%横坐标为Hz 
kx=round([1:M]*50/fmin);
%各次谐波对应的下标(从0开始) 
for m=1:M 
    km(m)=searchpeaks(amp,kx(m)+1);			%km为谱峰(从1开始) 
    if(amp(km(m)+1)<amp(km(m)-1)) 
        km(m)=km(m)-1; 
    end 
    beta(m)=amp(km(m)+1)./amp(km(m)); 
    delta(m)=(2*beta(m)-1)./(1+beta(m)); 
end 
fx=(km-1+delta)*fmin;							%估计频率 
hrx=amp(km)*2.*pi.*delta.*(1-delta.*delta)./sin(pi*delta);
degx=pha(km)-delta.*180/N*(N-1);				%估计相位 
degx=mod(degx,360);							%调整到0-360度 
efx=(fx-f)./f*100;								%频率误差 
ehr=(hrx-hr)./hr*100;							%幅度误差 
edeg=(degx-deg);								%相位误差 
% 结果输出: 
subplot(2,2,1);
%画出采样序列 
plot(t,x); 
hold on; 
plot(t,x.*w,'r');
%加窗波形 
hold off; 
xlabel('x(k)'); 
title('原信号和加窗信号 '); 
subplot(2,2,2);
%画出FFT分析结果 
stem(xaxis,amp,'.r'); 
xlabel('频率'); 
title('幅频结果'); 
subplot(2,2,4); 
stem(xaxis,pha,'.r'); 
xlabel('角频率'); 
title('相频结果'); 
subplot(2,2,3); 
stem(ehr); 
title('幅度误差(%)'); 
%文本输出 
fid=fopen('result.txt','w'); 
fprintf(fid,'原始数据:f1=%6.1fHz, N=%.f,  fs=%.f \r\n\r\n',f1,N,fs); 
fprintf(fid,'谐波次数      1      2      3      4      5      6     7\r\n'); 
fprintf(fid,'设定频率 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',f); 
fprintf(fid,'估计频率 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',fx); 
fprintf(fid,'误差(%%)  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n\r\n',efx); 
fprintf(fid,'设定幅值 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',hr); 
fprintf(fid,'估计幅值 %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n',hrx); 
fprintf(fid,'误差(%%)  %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f\r\n\r\n',ehr); 
fprintf(fid,'设定相位 %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n',deg); 
fprintf(fid,'估计相位 %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n',degx); 
fprintf(fid,'误差(度)  %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f\r\n\r\n',edeg); 
%其他数据 
fprintf(fid,'谱峰位置理论值:\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',[1:M]*f1/fmin); 
fprintf(fid,'谱峰位置估计值:\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',km-1+delta); 
fprintf(fid,'误差(%%)\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',((km-1+delta)-[1:M]*f1/fmin)./([1:M]*f1/fmin)*100); 
fprintf(fid,'delta     :\r\n %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f %6.4f\r\n',delta); 
fclose(fid);

%运行过程中的调用子程序:
function index1=searchpeaks(x,index) 
%在数组中寻找最大值对应的下标 
%x为数组,index 为给定的下标(index不能取最前或最后2个下标),在前后2个数中(共5个数)查找最大值和紧邻的次最大值 
% indexmax 返回两个谱峰位置中的前一个谱峰对应的下标 
index1=index-2; 
for k=-1:2 
    if(x(index+k)>x(index1)) 
        index1=index+k; 
    end 
end 
if x(index1-1)>x(index1+1) 
    index1=index1-1; 
end


#result.txt输出结果
原始数据:f1=  50.1Hz, N=2048,  fs=10000 

谐波次数      1      2      3      4      5      6     7
设定频率 50.100 100.200 150.300 200.400 250.500 300.600 350.700
估计频率 50.100 78.819 150.302 181.252 250.499 279.138 350.701
误差(%)  -0.000 -21.338  0.001 -9.555 -0.000 -7.140  0.000

设定幅值 35.355  0.000  2.482  0.000  1.252  0.000  1.125
估计幅值 35.356  0.046  2.482  0.002  1.252  0.002  1.125
误差(%)   0.001    Inf  0.009    Inf  0.004    Inf  0.003

设定相位  10.00   0.00  40.00   0.00  70.00   0.00 110.00
估计相位  10.03  31.67  39.97 338.35  70.06 329.86 110.05
误差()    0.03  31.67  -0.03 338.35   0.06 329.86   0.05

谱峰位置理论值:
 10.2605 20.5210 30.7814 41.0419 51.3024 61.5629 71.8234
谱峰位置估计值:
 10.2605 16.1421 30.7818 37.1203 51.3022 57.1675 71.8235
误差(%)
 -0.0002 -21.3385 0.0012 -9.5551 -0.0004 -7.1396 0.0002
delta     :
 0.2605 0.1421 0.7818 0.1203 0.3022 0.1675 0.8235

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-wL0WYF6F-1573836815517)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanning_han.bmp)]

汉明窗FIR滤波器设计:

汉宁窗(hanming window)又称改进升余弦窗,窗函数为:
wHn(n)=[0.540.46cos(2πnN1)]RN(n) w_{Hn}(n) = [0.54-0.46*cos(\frac{2\pi*n}{N-1})]*R_N(n)
幅值函数为:
WHn(w)=0.54WR(w)+0.23[WR(w2πN1)+0.23WR(w+2πN1)] W_{Hn}(w) = 0.54W_R(w) + 0.23[W_R(w-\frac{2\pi}{N-1})+0.23W_R(w+\frac{2\pi}{N-1})]
汉明窗主瓣宽度与汉宁窗相同,8π/N99.96%8{\pi}/N,99.96\%的能量集中在主瓣,第一瓣比主瓣低41dB。

在Matlab中,实现汉宁窗的函数为hanming ,其调用格式如下:

w=hanming(N)

案例分析1:设计一个汉明窗低通滤波器:

%语音信号设计一个汉明窗低通滤波器:
[x,FS,bits]=wavread('C:\Windows\Media\Windows Ringout');
x=x(:,1);
figure(1);
subplot(211);plot(x);
title('语音信号时域波形图')
xlabel('n');ylabel('h(n)')
grid on
y=fft(x,1000);
f=(FS/1000)*[1:1000];  
subplot(212);
plot(f(1:300),abs(y(1:300)));
title('语音信号频谱图');
xlabel('频率');ylabel('幅度')
grid on
%产生噪声信号并加到语音信号
t=0:length(x)-1;
zs0=0.05*cos(2*pi*10000*t/1024);
zs=[zeros(0,20000),zs0];
figure(2);
subplot(211)
plot(zs)
title('噪声信号波形');
xlabel('n');ylabel('h(n)')
grid on
zs1=fft(zs,1200);
subplot(212)
plot(f(1:600),abs(zs1(1:600)));
title('噪声信号频谱');
xlabel('频率');ylabel('幅度')
grid on
x1=x+zs';
%sound(x1,FS,bits);
y1=fft(x1,1200);
figure(3);
subplot(211);plot(x1);
title('加入噪声后的信号波形');
xlabel('n');ylabel('h(n)')
grid on
subplot(212);
plot(f(1:600),abs(y1(1:600)));
title('加入噪声后的信号频谱');
xlabel('频率');ylabel('幅度')
grid on
%滤波
fp=7500;
fc=8500; 
wp=2*pi*fp/FS;
ws=2*pi*fc/FS;
Bt=ws-wp; 
N0=ceil(6.2*pi/Bt);     
N=N0+mod(N0+1,2);
wc=(wp+ws)/2/pi;         
hn=fir1(N-1,wc,hamming(N)); 
X=conv(hn,x);           
X1=fft(X,1200);
figure(4);
subplot(211);
plot(X);
title('滤波后的信号波形');
xlabel('n');ylabel('h(n)')
grid on
subplot(212);
plot(f(1:600),abs(X1(1:600))); 
title('滤波后的信号频谱')
xlabel('频率');ylabel('幅度')
grid on

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-IRdanKny-1573836815518)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanming_lowp_1.bmp)]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zpDvFiHo-1573836815519)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanming_lowp_noise.bmp)]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-bBaOdODs-1573836815521)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanming_lowp_addnoise.bmp)]

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-iBooX8kL-1573836815523)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanming_lowp_hanfilter.bmp)]

案例分析2:已知连续信号为xa(t)=cos(100πt)+sin(100πt)+cos(50πt)x_a(t)=cos(100{\pi}t) +sin(100{\pi}t)+cos(50{\pi}t),用DFT分析其中xa(t)x_a(t)的频谱结构,选择不同的截取长度TpT_p。观察存在的截断效应,试用加窗的方法减少谱间干扰。

clear;close all
fs=400;T=1/fs;						%采样频率和采样间隔
Tp=0.04;N=Tp*fs;					%采样点数N
N1=[N,4*N,8*N];					%设定三种截取长度 
for m=1:3
    n=1:N1(m);
    xn=cos(100*pi*n*T)+ sin(200*pi*n*T)+ cos(50*pi*n*T);
    Xk=fft(xn,4096);
fk=[0:4095]/4096/T;
subplot(3,2,2*m-1);plot(fk,abs(Xk)/max(abs(Xk)));
if m==1 title('矩形窗截取');
end
end
%hamming窗截断
for m=1:3
    n=1:N1(m);
    wn=hamming(N1(m));
    xn=cos(200*pi*n*T)+ sin(100*pi*n*T)+ cos(50*pi*n*T).*wn';
    Xk=fft(xn,4096);
fk=[0:4095]/4096/T;
subplot(3,2,2*m);plot(fk,abs(Xk)/max(abs(Xk)));
if m==1 title('hamming窗截取');
end
end

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-39eMxIRF-1573836815524)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\hanming_dft.bmp)]

比较矩形窗与汉明窗的普分析结果可见,矩形窗比用汉明窗分辨率高(泄露小),但是谱间干扰大。因此汉明窗是以牺牲分辨率来换取谱间干扰降低。

布莱克曼FIR滤波器设计:

布莱克曼(Blackman window)的窗函数为:
wBl(n)=[0.420.5cos(2πnN1)+0.08cos(4πnN1)]RN(n) w_{Bl}(n) = [0.42-0.5*cos(\frac{2\pi*n}{N-1})+0.08*cos(\frac{4\pi*n}{N-1})]*R_N(n)
幅值函数为:
WHn(w)=0.42WR(w)+0.25[WR(w2πN1)+WR(w+2πN1)] W_{Hn}(w) = 0.42W_R(w) + 0.25[W_R(w-\frac{2\pi}{N-1})+W_R(w+\frac{2\pi}{N-1})]

+0.04[WR(w4πN1+WR(w+4πN1)] +0.04[W_R(w-\frac{4\pi}{N-1}+W_R(w+\frac{4\pi}{N-1})]

布莱克曼窗幅度函数由5部分相加而成,5部分相加的结果使得旁瓣得到进一步抵消,阻带衰减加大而过渡带加大到12π/N12{\pi}/N

在Matlab中,实现布莱克曼窗的函数为blackman ,其调用格式如下:

w=blackman(N);

:案例:用窗函数法设计数字带通滤波器。下阻带边缘:Ws1=0.3πAs=65dBW_{s1}=0.3{\pi},A_s=65dB,下通带边缘:Wp1=0.4πRp=1dBW_{p1}=0.4{\pi},R_p=1dB,上通带边缘:Wp2=0.6πRp=1dBW_{p2}=0.6{\pi},R_p=1dB,上阻带边缘:Ws2=0.7πRp=65dBW_{s2}=0.7{\pi},R_p=65dB。根据窗函数最小阻带衰减的特性,以及参照窗函数的基本参数表,选择布莱克曼窗可以达到75dB的最小阻带衰减,其过渡带为11π/N11\pi/N

clear all;
wp1=0.4*pi;
wp2=0.6*pi;
ws1=0.3*pi;
ws2=0.7*pi;
As=65;
tr_width=min((wp1-ws1),(ws2-wp2)); 					%过渡带宽度 
M=ceil(11*pi/tr_width)+1							%滤波器长度
n=[0:1:M-1];
wc1=(ws1+wp1)/2;									%理想带通滤波器的下截止频率
wc2=(ws2+wp2)/2;									%理想带通滤波器的上截止频率
hd=ideal_lp(wc2,M)-ideal_lp(wc1,M);
w_bla=(blackman(M))';								%布莱克曼窗
h=hd.*w_bla;
%截取得到实际的单位脉冲响应
[db,mag,pha,grd,w]=freqz_m(h,[1]);
%计算实际滤波器的幅度响应
delta_w=2*pi/1000;
Rp=-min(db(wp1/delta_w+1:1:wp2/delta_w))
%实际通带纹波
As=-round(max(db(ws2/delta_w+1:1:501)))
As=75
subplot(2,2,1);
stem(n,hd);
title('理想单位脉冲响应hd(n)')
axis([0 M-1 -0.4 0.5]);
xlabel('n');
ylabel('hd(n)')
grid on;
subplot(2,2,2);
stem(n,w_bla);
title('布莱克曼窗w(n)')
axis([0 M-1 0 1.1]);
xlabel('n');
ylabel('w(n)')
grid on;
subplot(2,2,3);
stem(n,h);
title('实际单位脉冲响应hd(n)')
axis([0 M-1 -0.4 0.5]);
xlabel('n');
ylabel('h(n)')
grid on;
subplot(2,2,4);
plot(w/pi,db);
axis([0 1 -150 10]);
title('幅度响应(dB)');
grid on;
xlabel('频率单位:pi');
ylabel('分贝数')

%调用小程序设计1:
function [db,mag,pha,grd,w] = freqz_m(b,a);
% Modified version of freqz subroutine
% ------------------------------------
% [db,mag,pha,grd,w] = freqz_m(b,a);
%  db = Relative magnitude in dB computed over 0 to pi radians
% mag = absolute magnitude computed over 0 to pi radians 
% pha = Phase response in radians over 0 to pi radians
% grd = Group delay over 0 to pi radians
%   w = 501 frequency samples between 0 to pi radians
%   b = numerator polynomial of H(z)   (for FIR: b=h)
%   a = denominator polynomial of H(z) (for FIR: a=[1])
%
[H,w] = freqz(b,a,1000,'whole');
    H = (H(1:1:501))'; w = (w(1:1:501))';
  mag = abs(H);
   db = 20*log10((mag+eps)/max(mag));
  pha = angle(H);
%  pha = unwrap(angle(H));
  grd = grpdelay(b,a,w);
%  grd = diff(pha);
%  grd = [grd(1) grd];
%  grd = [0 grd(1:1:500); grd; grd(2:1:501) 0];
%  grd = median(grd)*500/pi;

%调用小程序设计2:
function hd=ideal_lp(wc,M);
%计算理想低通滤波器的脉冲响应
%[hd]=ideal_lp(wc,M)
%hd=理想脉冲响应0到M-1
%wc=截止频率
% M=理想滤波器的长度
alpha=(M-1)/2;
n=[0:1:(M-1)];
m=n-alpha+eps;
%加上一个很小的值eps避免除以0的错误情况出现
hd=sin(wc*m)./(pi*m);


[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-52LsWSri-1573836815526)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\blackmen_bandp.bmp)]

凯塞窗FIR滤波器设计:

凯塞-贝塞窗(Kaiser-Basel window)的窗函数为:
wk(n)=I0(β)I0(α)RN(n) w_{k}(n) = \frac{I_0( β )}{I_0( \alpha )}*R_N(n)
式中,β=α(1(2nN1)1)2\beta= {\alpha}{\sqrt{(1-(\frac{2n}{N-1})-1)^2}}I0(x)I_0( x )是零阶第一类修正贝塞函数,可用下面级数计算:
I0(x)=1+k=1+1k!(x2)k2 I_0( x ) = 1+\sum_{k=1}^{+ ∞}(\frac{1}{k!}({\frac{x}{2}})^{k})^{2}
I0(x)I_0( x )q取15-25项就可以满足精度要求。通常α\alpha用以控制窗的形状,α\alpha加大,主瓣加宽,旁瓣减小,典型数据4<α<94<\alpha<9。当α=5.44\alpha=5.44s时,窗函数接近汉明窗;当α=7.865\alpha=7.865s时,窗函数接近于布莱克曼窗。其幅值函数为:
Wk(w)=wk(0)+2n=1(N1)2(wk(n)cos(wn)) W_{k}(w) =w_k(0) + 2\sum_{n=1}^{\frac{(N-1)}{2}}(w_k(n)cos(wn))
在Matlab中,实现汉宁窗的函数为kaiser,其调用格式如下:

w=kaiser(N);

在Matlab中设计标准响应FIR滤波器可使用fir1函数。fir1函数以经典方法实现加窗性相位FIR滤波器的设计,它可以设计出标准的低通、高通、带通、带阻滤波器。fir1函数用法为:

b = fir1(n,Wn,‘ftype’,wimdow)

各个参数的含义如下:

  • b -滤波器系数,n-滤波器阶数。
  • Wn -截止频率,0<=Wn<=10<=W_n<=1,Wn=1W_n=1对应于采样频率的一半。当设计带通和带阻滤波器时,Wn=[W1,W2],W1<w<W2W_n =[W_1,W_2],W_1<w<W_2
  • ftype -当指定ftype时,可设计高通和带阻滤波器。ftype=hight时,设计高通FIR滤波器;ftype=stop时设计带阻FIR滤波器。低通和带通FIR滤波器无需输入ftype参数。
  • window–窗函数。窗函数的长度应等于FIR滤波器系数的个数,即阶数n+1。

案例分析:利用凯塞窗函数设计一个带通滤波器,上截止频率2500Hz,下截止频率1000Hz,过渡带宽200Hz,通带纹波允许差0.1,带阻纹波不大于允差0.02dB,通带幅值为1。

Fs=8000;N=216;
fcuts=[1000 1200 2300 2500];
mags=[0 1 0];
devs=[0.02 0.1 0.02];
[n,Wn,beta,ftype]=kaiserord(fcuts,mags,devs,Fs);
n=n+rem(n,2);
hh=fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
[H,f]=freqz(hh,1,N,Fs);
plot(f,abs(H));
xlabel('频率 (Hz)');
ylabel('幅值|H(f)|');
grid on;

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zSlKyrYh-1573836815527)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\kaser_bandp.bmp)]

窗函数设计法:

根据前面几节的分析:设计一个FIR低通滤波器通常按照下面的步骤执行:

  1. 根据滤波器设计要求,确定滤波器的过渡带宽和阻带衰减要求,选择合适窗函数的类型并进行估计窗函数的宽度N。
  2. 根据所求的理想滤波器求出单位脉冲响应hd(n)h_d(n)
  3. 根据求得的h(n)求出其频率响应。
  4. 根据频率响应验证是否满足技术指标。
  5. 若不满足指标要求,则应调整窗函数类型或者长度,然后重复(1),(2),(3)(4)步,直到满足要求为止。

注意:matlab中数据通常是以列向量形式存储的,所以两个向量相乘必须进行转置。计算滤波器的单位脉冲响应h(n),根据窗函数设计理论h(n)=hd(n)w(n)h(n)=h_d(n)*w(n),在matlab中用语句hn=hd*wd实现h(n).

窗函数设计法程序设计如下:

function [h]=usefir1(mode,n,fp,fs,window,r,sample)
% mode:模式(1--高通; 2--低通; 3--带通; 4--带阻)
% n:阶数, 加窗的点数为阶数加1
% fp:高通和低通时指示截止频率, 带通和带阻时指示下限频率
% fs:带通和带阻时指示上限频率
% window:加窗(1--矩形窗; 2--三角窗; 3--巴特窗; 4--汉明窗; 
%5--汉宁窗; 6--布莱克曼窗; 7--凯泽窗; 8--契比雪夫窗)
% r代表加chebyshev窗的r值和加kaiser窗时的beta值
% sample:采样率
% h:返回设计好的FIR滤波器系数
if window==1 w=boxcar(n+1);
end
if window==2 w=triang(n+1);end
if window==3 w=bartlett(n+1);end
if window==4 w=hamming(n+1);end
if window==5 w=hanning(n+1);end
if window==6 w=blackman(n+1);end
if window==7 w=kaiser(n+1,r);end
if window==8 w=chebwin(n+1,r);
end
wp=2*fp/sample;
ws=2*fs/sample;
if mode==1 h=fir1(n,wp,'high',w);
end
if mode==2 h=fir1(n,wp,'low',w);
end
if mode==3 h=fir1(n,[wp,ws],w);
end
if mode==4 h=fir1(n,[wp,ws],'stop',w);
end
m=0:n;
subplot(131);
plot(m,h);grid on;	
title('冲激响应');
axis([0 n 1.1*min(h) 1.1*max(h)]);
ylabel('h(n)');xlabel('n');
freq_response=freqz(h,1);
magnitude=20*log10(abs(freq_response));
m=0:511; f=m*sample/(2*511);
subplot(132);
plot(f,magnitude);grid on;
title('幅频特性');
axis([0 sample/2 1.1*min(magnitude) 1.1*max(magnitude)]);
ylabel('f幅值');xlabel('频率');
phase=angle(freq_response);
subplot(133);plot(f,phase);grid on;
title('相频特性');
axis([0 sample/2 1.1*min(phase) 1.1*max(phase)]);
ylabel('相位');xlabel('频率');

案例分析:假设需要设计一个40阶的带通FIR滤波器,采用汉明窗,采样频率为10kHz,两个截止频率分别为2kHz和3kHz,则需要在Matlab的命令行窗口输入:

h=usefir1(3,60,2000,3000,4,2,10000);

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-8Mmw8ohg-1573836815529)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\usefir1.bmp)]

2019-11-14 23:34:56 qq_17119267 阅读数 321

[Matlab]FIR滤波器设计:(FIR滤波器的结构)

FIR(Finite Impulse Response)滤波器:有限长单位冲激响应滤波器,又称为非递归型滤波器,是一种在数字信号领域应用非常广泛的基础型滤波器。FIR滤波器的特点是能够在输入具体任意幅频特性的数字信号后,保证输出数字信号的相频特性仍然保持严格线性。另外FIR滤波器是数字信号处理系统中最基本的元件,它可以在保证任意幅频特性的同时具有严格的线性相频特性,同时具有有限长的脉冲采样响应特性,因而滤波器是稳定的系统。因此,FIR滤波器的应用要远胜于IIR滤波器,在通信、图像处理、模式识别等领域都有着广泛而举足轻重的作用。

FIR滤波器的结构:

​ 尽管IIR数字滤波器的设计理论已经非常成熟、经典,但是IIR滤波器有一个自身的中要缺陷:相位特性通常是非线性的。

  • IIR滤波器的非线性原理:在IIR滤波器的滤波器设计中,只是对滤波器的幅频特性进行了研究,并获得了良好的幅频特性。但对相频特性却没有考虑,所以IIR滤波器的相频特性通常是非线性的。
  • IIR滤波器的应用问题:由于IIR滤波器的相位特性通常是非线性的,所以它适合用于在对系统相位特性要求不严格的场合。但是另一方面,由于一个具有线性相位特性的滤波器可以保证在滤波器通带内信号传输不失真,所以在许多领域需要滤波器具有严格得线性相位特性,比如图像处理以及数据传输等。但是如果要使用IIR滤波器实现线性的相位特性,则必须对其相位特性用全通滤波器进行校正,其结果使得滤波器设计变得非常复杂,实现困难成本提高。所以,要实现具有线性相位的数字滤波器,必须另辟蹊径。
  • 有限脉冲响应系统的单位脉冲响应h(s)h(s)为有限长序列,系统函数H(z)H(z)在有限z平面上不存在几点,其运算结构中没有反馈支路,即没有环路。所以,有限脉冲响应滤波器可以设计成在整个频率范围内均可以提供精确的线性相位,而且总是可以独立于滤波器系数保持有限长输入有限长输出的稳定。因此,这样的滤波器在很多领域是首选的。
FIR滤波器的特点

有限长单位冲激响应(FIR)滤波器有以下特点:

  1. 系统的单位冲激响应h (n)在有限个n值处不为零
  2. 系统函数H(z)在|z|>0处收敛,极点全部在z = 0处(因果系统
  3. 结构上主要是非递归结构,没有输出到输入的反馈
FIR滤波器的基本结构:
  • 直接型:

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-xTYGXjLS-1573745391892)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\截图20191114230732.png)]

  • 级联型:

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-vfmwWIuD-1573745391896)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\截图20191114230814.png)]

特点:

  1. 每个基本节控制一对零点便于控制滤波器的传输零点
  2. 系数比直接型多所需的乘法运算多。

直接型与级联型:

clear all;
n=0:10;
N=30;
b=0.9.^n;
delta=impseq(0,0,N);%脉冲系列生成器
h=filter(b,1,delta);%其中b是分子系数向量,a是分母系数向量。
x=[ones(1,5),zeros(1,N-5)];
y=filter(b,1,x);
subplot(2,2,1);stem(h);%绘制火柴图
title('直接型h(n)');
subplot(2,2,2);stem(y);
title('直接型y(n)');
[b0,B,A]=dir2cas(b,1);
h=casfilter(b0,B,A,delta);%级联型单位脉冲响应
y=casfilter(b0,B,A,x);
subplot(2,2,3);stem(h);
title('级联型h(n)');
subplot(2,2,4);stem(y);
title('级联型y(n)');

%%impseq函数设计生成m文件放入当前文件夹:
function [x,n] = impseq(n0,n1,n2)
% Generates x(n) = delta(n-n0); n1 <= n,n0 <= n2
% ----------------------------------------------
% [x,n] = impseq(n0,n1,n2)
%
if ((n0 < n1) | (n0 > n2) | (n1 > n2))
	error('arguments must satisfy n1 <= n0 <= n2')
end
n = [n1:n2];
%x = [zeros(1,(n0-n1)), 1, zeros(1,(n2-n0))];
x = [(n-n0) == 0];

%%dir2cas函数设计生成m文件放入当前文件夹:
function [b0,B,A] = dir2cas(b,a)
% 直接型到级联型的型式转换(复数对型)
% [b0,B,A] = dir2cas(b,a)
% b = 直接型的分子多项式系数
% a = 直接型的分母多项式系数
% b0 = 增益系数
% B = 包含各bk的K乘3维实系数矩阵
% A = 包含各ak的K乘3维实系数矩阵
% compute gain coefficient b0
b0 = b(1); b = b/b0;
a0 = a(1); a = a/a0;
b0 = b0/a0;
%
M = length(b); N = length(a);
if N > M
	b = [b zeros(1,N-M)];
elseif M > N
	a = [a zeros(1,M-N)]; N = M;
else
	NM = 0;
end
%
K = floor(N/2); B = zeros(K,3); A = zeros(K,3);
if K*2 == N;
	b = [b 0];
	a = [a 0];
end
%        
broots = cplxpair(roots(b));
aroots = cplxpair(roots(a));
for i=1:2:2*K
	Brow = broots(i:1:i+1,:);
	Brow = real(poly(Brow));
	B(fix((i+1)/2),:) = Brow;
	Arow = aroots(i:1:i+1,:);
	Arow = real(poly(Arow));
	A(fix((i+1)/2),:) = Arow;
end     

%%casfilter函数设计生成m文件放入当前文件夹:
function y = casfilter(b0,B,A,x)
% CASCADE form realization of IIR and FIR filters
% -----------------------------------------------
% y = casfiltr(b0,B,A,x);
%  y = output sequence
% b0 = gain coefficient of CASCADE form
%  B = K by 3 matrix of real coefficients containing bk's
%  A = K by 3 matrix of real coefficients containing ak's
%  x = input sequence
%
[K,L] = size(B);
N = length(x);
w = zeros(K+1,N);
w(1,:) = x;
for i = 1:1:K
        w(i+1,:) = filter(B(i,:),A(i,:),w(i,:));
end
y = b0*w(K+1,:);

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-wyQ1ZD5y-1573745391900)(G:\研究生\项目小组任务\程序设计\滤波器设计\first.bmp)]

  • 快速卷积型:

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-Ss7vEm1J-1573745391897)(G:\研究生\项目小组任务\笔记\第四周和第五周笔记\截图20191114230922.png)]

运用FFT实现有限长序列x(n)x(n)h(x)h(x)的线性卷积,则可以得到FIR滤波器的快速卷积结构如图所示,对于x(n)x(n)的无限长序列的一般情况,可以采用重叠相加法或者重叠保留法实现FIR滤波器的快速卷积结构

  • 频率采样型:

FIR系统直接型转换为频率采样型结构的Matlab实现源代码:

function [C,B,A] = dir2fs(h)
% 直接型到频率采样型的转换
% [C,B,A] = dir2fs(h)
% C = 包含各并行部分增益的行向量
% B = 包含按行排列的分子系数矩阵
% A = 包含按行排列的分母系数矩阵
% h =  FIR滤波器的脉冲响应向量
M = length(h);
H = fft(h,M);
magH = abs(H); phaH = angle(H)';
% check even or odd M
if (M == 2*floor(M/2))
      L = M/2-1;   %  M为偶数 
     A1 = [1,-1,0;1,1,0];
     C1 = [real(H(1)),real(H(L+2))];
 else
      L = (M-1)/2; % M is odd
     A1 = [1,-1,0];
     C1 = [real(H(1))];
 end
k = [1:L]';
% 初始化 B 和 A 数组
B = zeros(L,2); A = ones(L,3);
% 计算分母系数
A(1:L,2) = -2*cos(2*pi*k/M); A = [A;A1];
% 计算分子系数
B(1:L,1) = cos(phaH(2:L+1));
B(1:L,2) = -cos(phaH(2:L+1)-(2*pi*k/M));
% 计算增益系数
C = [2*magH(2:L+1),C1]';

案例分析:

%[例6-1]
%%利用频率采样法设计一个低通FIR数字低通滤波器,其理想频率特性是矩形的
%pi=3.14
%给定抽样频率为Ωs=2*pi*1.5*10^4(rad/sec),通带截至频率为Ωp=2*pi*1.6*10^3(rad/sec),
%阻带截至频率为Ωst=2*pi*3.1*10^3(rad/sec),通带波动<=1dB ,阻带衰减>=50dB。
%通带的截至频率为:w_p = Ωp/f_s = 2*pi Ωp/Ωs = 0.213*pi
%阻带的起始频率为:w_st = Ωst/f_s = 2*pi Ωfs/Ωs = 0.413*pi
%理想低通截至频率:Ωc = 1/2(Ωp+Ωst) = 2*pi*2.35*10^3(rad/sec)
%其对应的数字频率:w_c = 2*pi Ωc/Ωs 0.313*pi
%过渡带带宽为:△w = w_st - w_p =0.2*pi  由于△w = (2*pi/N) *3 所以抽样频率为30rad/sec
% 直接型到频率采样型的转换
close all;
clear;
N=30;
H=[ones(1,4),zeros(1,22),ones(1,4)];
H(1,5)=0.5886;H(1,26)=0.5886;H(1,6)=0.1065;H(1,25)=0.1065;
k=0:(N/2-1);k1=(N/2+1):(N-1);k2=0;
A=[exp(-j*pi*k*(N-1)/N),exp(-j*pi*k2*(N-1)/N),exp(j*pi*(N-k1)*(N-1)/N)];
HK=H.*A;
h=ifft(HK);
fs=15000;
[c,f3]=freqz(h,1);
f3=f3/pi*fs/2;
subplot(221);
plot(f3,20*log10(abs(c)));
title('频谱特性');
xlabel('频率/HZ');ylabel('衰减/dB');
grid on;
subplot(222);
title('输入采样波形');
stem(real(h),'.');
line([0,35],[0,0]);
xlabel('n');ylabel('Real(h(n))');
grid on;
t=(0:100)/fs;
W=sin(2*pi*t*750)+sin(2*pi*t*3000)+sin(2*pi*t*6500);
q=filter(h,1,W);
[a,f1]=freqz(W);
f1=f1/pi*fs/2;
[b,f2]=freqz(q);
f2=f2/pi*fs/2;
subplot(223);
plot(f1,abs(a));
title('输入波形频谱图');
xlabel('频率');ylabel('幅度')
grid on;
subplot(224);
plot(f2,abs(b));
title('输出波形频谱图');
xlabel('频率');ylabel('幅度')
grid on;

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-03u5QQ1I-1573745391902)(G:\研究生\项目小组任务\程序设计\滤波器设计\secord.bmp)]

2019-06-16 13:06:49 BIGFatming 阅读数 288

基于C#的窗函数法低通FIR滤波器

摘要

基于Visual Studio 2015 开发环境,使用C#编程语言,运用窗函数法构造了FIR低通滤波器。并通过MATLAB对其滤波效果进行了试验。
由于凯泽窗的数学模型过于复杂,笔者能力有限,故无法在本文中给出凯泽窗的窗函数代码。

源代码
1、矩形窗函数

namespace Low_Pass_FIR_Fliter
{
    class Boxcar
    {
        public int N = 0;
        public Boxcar (double Wp, double Ws)
        {
            int i;
            double n = (0.9 * 2 * Math.PI) / (Ws - Wp);//计算滤波器长度
            for (i = 0; i < n; i++) ;
            N = i;//得到滤波器长度
        }
        public double[] Get_Win()//返回N点矩形窗序列
        {
            int n;
            double[] wd = new double[N];
            for (n = 0; n < N; n++)
            {
                wd[n] = 1;
            }
            return wd;
        }
    }
}

2、三角窗函数

namespace Low_Pass_FIR_Fliter
{
    class Bartlett
    {
        public int N = 0;
        public Bartlett (double Wp, double Ws)
        {
            int i;
            double n = (2.1 * 2 * Math.PI) / (Ws - Wp);//计算滤波器长度
            for (i = 0; i < n; i++) ;
            N = i;//得到滤波器长度
        }
        public double[] Get_Win()//返回N点三角窗序列
        {
            int n;
            double[] wd = new double[N];
            for (n = 0; n < N; n++)
            {
                if(n<=((N-1)/2))
                {
                    wd[n] = (2 * n / (N - 1));
                }
                if(n > ((N - 1) / 2))
                {
                    wd[n] = 2 - (2 * n / (N - 1));
                }
            }
            return wd;
        }
    }

3、汉宁窗

namespace Low_Pass_FIR_Fliter
{
    class Hanning
    {
        public int N = 0;
        public Hanning(double Wp, double Ws)
        {
            int i;
            double n = (3.1 * 2 * Math.PI) / (Ws - Wp);//计算滤波器长度
            for (i = 0; i < n; i++) ;
            N = i;//得到滤波器长度
        }
        public double[] Get_Win()//返回N点汉宁窗序列
        {
            int n;
            double[] wd = new double[N];
            for (n = 0; n < N; n++)
            {
                double b = Math.Cos((2 * Math.PI*n) / (N - 1));
                double res = (0.5 - 0.5 * b);
                wd[n] = res;
            }
            return wd;
        }
    }
}

4、汉明窗

namespace Low_Pass_FIR_Fliter
{
    public class Hamming
    {
        public int N = 0;
        public Hamming(double Wp,double Ws)
        {
            int i;
            double n = (3.3 * 2 * Math.PI) / (Ws-Wp);//计算滤波器长度
            for (i = 0; i < n; i++) ;
            N = i ;//得到滤波器长度
        }
        public double [] Get_Win()//返回N点汉明窗序列
        {
            int n ;
            double[] wd=new double[N];
            for(n=0;n<N;n++)
            {
                double b = Math.Cos((2*n * Math.PI) / (N-1));
                double res = (0.54 - 0.46 * b);
                wd[n] = res;
            }
            return wd;
        }
    }
}

5、布莱克曼窗

namespace Low_Pass_FIR_Fliter
{
    class Blackman
    {
        public int N = 0;
        public Blackman (double Wp, double Ws)
        {
            int i;
            double n = (5.5 * 2 * Math.PI) / (Ws - Wp);//计算滤波器长度
            for (i = 0; i < n; i++) ;
            N = i;//得到滤波器长度
        }
        public double[] Get_Win()//返回N点布莱克曼窗序列
        {
            int n;
            double[] wd = new double[N];
            for (n = 0; n < N; n++)
            {
                double a = Math.Cos((2 * n * Math.PI) / (N - 1));
                double b = Math.Cos((4 * n * Math.PI) / (N - 1));
                double res = (0.42 - 0.5 * a + 0.08 * b);
                wd[n] = res;
            }
            return wd;
        }
    }
}

6、单位冲激响应

namespace Low_Pass_FIR_Fliter
{
    class UnitImpulseReact
    {
        private double Wc;
        private int alpha;
        private int N;
        public UnitImpulseReact(double Wp,double Ws,int N)//构造单位冲激响应
        {
            this.Wc = 0.5 * (Wp + Ws);//取近似截止频率
            if (N % 2 != 0)
            {
                this.alpha = (N - 1) / 2;//计算相移常数
            }
            else
            {
                this.alpha = N / 2;
            }
            this.N = N;
        }
        public double[] Get_UIR()//获取单位冲激响应序列
        {
            double[] hd = new double[N];
            for(int n=0;n<N;n++)
            {
                double numerator = Math.Sin(Wc * (n - alpha));
                double denominator = Math.PI * (n - alpha);
                if (n == alpha)//当n=α时,取极限
                {
                    hd[n] = Wc / Math.PI;
                }
                else
                {
                    hd[n] = numerator / denominator;
                }
            }
            return hd;
        }
    }
}

使用MATLAB做正确性检验

大致思路:自定义一个15点序列,通过频谱图确定技术指标。分别通过MATLAB标准做法和自编程序对其进行滤波,若滤波结果大致一致,则可认为这个滤波器程序是正确的。

自定义序列 x[n]={ 16 12 2 -4 -1 1 -3 11 -27 18 15 -9 11 21 3} 该序列对应的频谱图为:
自定义序列的频谱特性
选取通带边界频率ωp=0.2π,阻带边界频率ωs=0.3π,阻带衰减>30dB。由此,选择汉宁窗。通过MATLAB程序构成的低通滤波器b1的频谱图为:
在这里插入图片描述
滤波后的序列y的频率特性为(将>0dB的部分局部放大)
在这里插入图片描述
同样的,在程序中基于预定的滤波器特性构建汉宁窗,并结合单位冲激响应得到滤波器的单位样值响应

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.IO;

namespace Low_Pass_FIR_Fliter
{
    class Program
    {
        static void Main(string[] args)
        {
            Hanning win = new Hanning(0.2 * Math.PI, 0.3 * Math.PI);//构造汉宁窗
            UnitImpulseReact UIR = new UnitImpulseReact(0.2 * Math.PI, 0.3 * Math.PI, win.N);//构造单位冲激响应
            double[] hd = new double[win.N];
            double[] w = new double[win.N];
            double[] h = new double[win.N];
            string[] res = new string[win.N];//定义数组
            w = win.Get_Win();//获取窗函数序列
            hd = UIR.Get_UIR();//获取单位冲激响应序列
            for(int i=0;i<win.N;i++)
            {
                h[i] = w[i] * hd[i];//得到滤波器序列
            }
            for (int i = 0; i < win.N; i++)
            {
                res[i] = Convert.ToString(h[i]);//将滤波器序列点数转化成字符串形式,准备写入文件
            }
            File.WriteAllLines(@"C:\VS2015\C#Programe\Low-Pass_FIR_Fliter\Low-Pass_FIR_Fliter\sequence.txt", res);//写入文件,用于MATLAB的后续处理
        }
    }
}

通过程序获得的滤波器b2的频率特性为:
在这里插入图片描述
程序滤波器滤波后的序列y1的频率特性为(同样将>0dB的部分局部放大)
在这里插入图片描述

通过波形比较可以认为程序生成的滤波器是成功的。

2019-04-16 21:14:02 weixin_38197667 阅读数 480

相关文章:
[Verilog上机实验题目1:8位数字显示的简易频率计]
[Verilog上机实验题目2:11位巴克码序列峰值检测器]
[Verilog上机实验题目3:FIR滤波器]
[Verilog上机实验题目4:哈夫曼编码器]

FIR滤波器的设计

 滤波器就是对特定的频率或者特定频率以外的频率进行消除的电路,被广泛应用于通信系统和信号处理系统中。从功能角度,数字滤波器对输入离散信号的数字代码进行运算处理,以达到滤除频带外信号的目的。

有限冲击响应(FIR)滤波器就是一种常用的数字滤波器,采用对已输入样值的加权形成它的输出。对于输入序列X[n]的FIR滤波器可用下图所示的结构示意图来表示,其中X[n]是输入数据流。各级的输入连接和输出连接被称为抽头,系数(b0,b1...)被称为抽头系数。一个M阶的FIR滤波器将会有M+1个抽头。通过移位寄存器用每个时钟边沿处 的数据流采样值乘以抽头系数,并将它们加起来形成输出。

FIR滤波器电路有两个主要功能模块:移位寄存器组模块Shift_reg用于储存串行进入滤波器的数据;乘加计算模块Caculater用于进行FIR计算。

顶层代码模块

module FIR(
	input clk,rst,
	input [3:0] Data_in,
	output[9:0] Data_out
 );
// 声明抽头类型为wire
	wire[3:0] 	samples_0, 	samples_1, 	samples_2, 	samples_3, 	samples_4, 	samples_5, 	samples_6, 	samples_7, 	samples_8; 
// 移位寄存器组模块实例化
	shift_reg U1(.Data_in(Data_in),.clk(clk),.rst(rst), 	.samples_0(samples_0), 	.samples_1(samples_1), 	.samples_2(samples_2), 	.samples_3(samples_3), 	.samples_4(samples_4), 	.samples_5(samples_5), 	.samples_6(samples_6), 	.samples_7(samples_7), 	.samples_8(samples_8)); 
// 乘加计算模块实例化
	caculator U2(.Data_out(Data_out), 	.samples_0(samples_0), 	.samples_1(samples_1), 	.samples_2(samples_2), 	.samples_3(samples_3), 	.samples_4(samples_4), 	.samples_5(samples_5), 	.samples_6(samples_6), 	.samples_7(samples_7), 	.samples_8(samples_8));  
endmodule  

shift_reg模块用于存储输入的数据流,本例中存储8个4位宽输入数据信号,也作为caculator模块的输入。

移位寄存器shift_reg模块

 module shift_reg( 	 input clk,rst, 	 input[3:0] Data_in, 	 output reg[3:0] 	 samples_0, 	samples_1, 	samples_2, 	samples_3, 	samples_4, 	samples_5, 	 samples_6, 	samples_7, 	samples_8  );  
 // 初始化清零,将数据一位一位移进,进行计算
 always@(posedge clk or negedge rst)
	 begin
		 if(!rst)
			 begin
				samples_0<=4'b0;
				samples_1<=4'b0;
				samples_2<=4'b0;
				samples_3<=4'b0;
				samples_4<=4'b0;
				samples_5<=4'b0;
				samples_6<=4'b0;
				samples_7<=4'b0;
				samples_8<=4'b0;
			 end
		 else 
			begin
				samples_0<=Data_in;
				samples_1<=samples_0;
				samples_2<=samples_1;
				samples_3<=samples_2;
				samples_4<=samples_3;
				samples_5<=samples_4;
				samples_6<=samples_5;
				samples_7<=samples_6;
				samples_8<=samples_7;
			end
	 end
 endmodule

caculator模块用于进行8输入信号与抽头系数的乘法和累加,并产生滤波后的信号data_out,抽头系数具有对称性,本例直接给出。

caculator模块

 module caculator(
	 input [3:0] 	samples_0, 	samples_1, 	samples_2, 	samples_3, 	samples_4, 	samples_5, 	samples_6, 	samples_7, 	samples_8, 	 
	 output [9:0] Data_out
);
	 wire[3:0] 	 out_t_1, 	 out_t_2, 	 out_t_3, 	 out_t_4;
 	 wire[7:0] 	 out1, 	 out2, 	 out3, 	 out4, 	 out5;
 // 将相同系数的2个输入信号先相加赋值给out_t
	 assign out_t_1 = samples_0 + samples_8;
	 assign out_t_2 = samples_1 + samples_7;
	 assign out_t_3 = samples_2 + samples_6;
	 assign out_t_4 = samples_3 + samples_5;
// 定义抽头系数为参数COEF
	parameter COEF1 = 4'b0010;
	parameter COEF2 = 4'b0011;
	parameter COEF3 = 4'b0110;
	parameter COEF4 = 4'b1010;
	parameter COEF5 = 4'b1100;
// 加法树乘法器实例   
	 mul_addtree U1(.mul_a(COEF1),.mul_b(out_t_1),.mul_out(out1));
	 mul_addtree U2(.mul_a(COEF2),.mul_b(out_t_2),.mul_out(out2));
	 mul_addtree U3(.mul_a(COEF3),.mul_b(out_t_3),.mul_out(out3));
	 mul_addtree U4(.mul_a(COEF4),.mul_b(out_t_4),.mul_out(out4));
	 mul_addtree U5(.mul_a(COEF5),.mul_b(samples_4),.mul_out(out5));
// 将所有的out加起来组成输出信号
	assign Data_out = out1 +   out2 + out3 + out4 + out5;
endmodule

 加法树乘法器模块

// 加法树乘法器模块   将out_t与对应的抽头系数相乘赋值给out
module mul_addtree(
	input [3:0]mul_a,mul_b,
	output [7:0] mul_out
);
	wire[7:0]stored0, stored1,stored2,stored3;
	wire[7:0] add01,add23;
	
	assign stored0 = mul_b[0] ? {4'b0,mul_a} : 8'b0;
	assign stored1 = mul_b[1] ? {3'b0,mul_a,1'b0} : 8'b0;
	assign stored2 = mul_b[2] ? {2'b0,mul_a,2'b0} : 8'b0;
	assign stored3 = mul_b[3] ? {1'b0,mul_a,3'b0} : 8'b0;
	
	assign add01 = stored0 + stored1;
	assign add23 = stored2 + stored3;
// 分成2级的原因是为了加快硬件运行速度
	assign mul_out = add01 + add23;
endmodule

测试代码如下

// 测试代码  16MHz周期是1/16M=0.0625us微秒=62.5ns纳秒,测试时间基准 1ns/10ps
 `timescale 1ns/10ps;
module 	FIR_tb;
	reg clk,rst;
	reg[3:0] Data_in;
	wire[9:0] Data_out;
	initial
		begin
			clk = 'b0;
			Data_in = 'b0;
			forever
				begin
					#31.25 clk <= ~clk ;	//时钟周期62.5ns
					#31.25 Data_in <= Data_in + 1'b1 ;	//斜坡波形
				end
		end
	
	initial
		begin
			rst = 1'b0;
			#150 rst = 1'b1;	//异步复位信号,一般等待2到3个周期再复位
		end
	
// FIR滤波器实例
FIR v1(.clk(clk),.rst(rst),.Data_in(Data_in),.Data_out(Data_out));
endmodule

仿真结果如下图所示,输入信号为斜坡信号,输出信号为经过FIR滤波器后的信号。

 

FIR滤波器

阅读数 1272

1---FIR滤波器简介

阅读数 20650

没有更多推荐了,返回首页