我借鉴了前辈的一段程序!是把音频信号加噪声然后用低通滤波器滤波!然后比较!
用的是Butterworth滤波器和blackman滤波器~但两个的频谱图错了~
我不知道是参数问题还是什么问题~大家指教下啊!
是不是要把音频的一些参数弄出来然后计算什么的啊!
我是一菜鸟~被老师一问就什么都不懂了~大家教教我啊!
blackman滤波器设计
Window=blackman(55); %长度为35的blackman窗Window
b=fir1(54,0.1,Window);%产生低通blackman滤波器
H=freqz(b,1,512);
figure(6);
plot(abs(H));
title('blackman滤波器的频率响应');
xlabel('频率/Hz');
ylabel('声强/db');
Butterworth滤波器设计
FS=1;
%通带、阻带截止频率
Fl=0.08;Fh=0.15;
%频率预畸
wp=(Fl/FS)*2*pi; %临界频率采用角频率表示
ws=(Fh/FS)*2*pi; %临界频率采用角频率表示
OmegaP=2*FS*tan(wp/2);
OmegaS=2*FS*tan(ws/2);
[k,Wn]=buttord(OmegaP,OmegaS,1.1,60,'s');
[b,a]=butter(k,Wn,'s');
%freqs(b,a) %设计模拟的
[bz,az]=bilinear(b,a,FS); %映射为数字的
% 绘制结果
H=freqz(bz,az,1024,FS,'whole');
figure(8);
plot(abs(H));
title('Butterworth滤波器的频率响应');
xlabel('频率/Hz');
ylabel('声强/db');
--------------------------
用的是Butterworth滤波器和blackman滤波器~但两个的频谱图错了~
我不知道是参数问题还是什么问题~大家指教下啊!
是不是要把音频的一些参数弄出来然后计算什么的啊!
我是一菜鸟~被老师一问就什么都不懂了~大家教教我啊!
blackman滤波器设计
Window=blackman(55); %长度为35的blackman窗Window
b=fir1(54,0.1,Window);%产生低通blackman滤波器
H=freqz(b,1,512);
figure(6);
plot(abs(H));
title('blackman滤波器的频率响应');
xlabel('频率/Hz');
ylabel('声强/db');
Butterworth滤波器设计
FS=1;
%通带、阻带截止频率
Fl=0.08;Fh=0.15;
%频率预畸
wp=(Fl/FS)*2*pi; %临界频率采用角频率表示
ws=(Fh/FS)*2*pi; %临界频率采用角频率表示
OmegaP=2*FS*tan(wp/2);
OmegaS=2*FS*tan(ws/2);
[k,Wn]=buttord(OmegaP,OmegaS,1.1,60,'s');
[b,a]=butter(k,Wn,'s');
%freqs(b,a) %设计模拟的
[bz,az]=bilinear(b,a,FS); %映射为数字的
% 绘制结果
H=freqz(bz,az,1024,FS,'whole');
figure(8);
plot(abs(H));
title('Butterworth滤波器的频率响应');
xlabel('频率/Hz');
ylabel('声强/db');
--------------------------
songzy41
blackman滤波器设计改为：
%blackman滤波器设计
Window=blackman(55); %长度为35的blackman窗Window
b=fir1(54,0.1,Window);%产生低通blackman滤波器
%freqz(b,1,512);
[H,w]=freqz(b,1,512);
plot(w/pi,abs(H)); grid;
title('blackman滤波器的频率响应');
xlabel('归一化频率');
ylabel('幅值');
给出的频率是归一化频率，Y轴不是声强/db，而是线性幅值。Butterworth滤波器设计改为
%Butterworth滤波器设计
FS=1;
T=1/FS;
%通带、阻带截止频率
Fl=0.08;Fh=0.15;
%频率预畸
wp=(Fl/FS)*2*pi; %临界频率采用角频率表示
ws=(Fh/FS)*2*pi; %临界频率采用角频率表示
%为双线性z变换进行频率变换
OmegaP=2*FS*tan(wp/2);
OmegaS=2*FS*tan(ws/2);
[k,Wn]=buttord(OmegaP,OmegaS,1,60,'s');
Wn=Wn/pi; %转换成数字滤波器的Wn
[b,a]=butter(k,Wn);
%freqs(b,a) %设计模拟的
%[bz,az]=bilinear(b,a,T); %映射为数字的
% 绘制结果
[H,f]=freqz(b,a,1024,'whole',FS);
%figure(8);
plot(f,abs(H)); grid;
title('Butterworth滤波器的频率响应');
xlabel('频率(Hz)');
ylabel('幅值');
figure(2);
[db, mag, pha, grd,w]=freqz_m(b,a);
plot(w*FS/2/pi,db); grid;
title('Butterworth滤波器的频率响应');
xlabel('频率(Hz)');
ylabel('幅值(dB)');
响应曲线如下
%blackman滤波器设计
Window=blackman(55); %长度为35的blackman窗Window
b=fir1(54,0.1,Window);%产生低通blackman滤波器
%freqz(b,1,512);
[H,w]=freqz(b,1,512);
plot(w/pi,abs(H)); grid;
title('blackman滤波器的频率响应');
xlabel('归一化频率');
ylabel('幅值');
给出的频率是归一化频率，Y轴不是声强/db，而是线性幅值。Butterworth滤波器设计改为
%Butterworth滤波器设计
FS=1;
T=1/FS;
%通带、阻带截止频率
Fl=0.08;Fh=0.15;
%频率预畸
wp=(Fl/FS)*2*pi; %临界频率采用角频率表示
ws=(Fh/FS)*2*pi; %临界频率采用角频率表示
%为双线性z变换进行频率变换
OmegaP=2*FS*tan(wp/2);
OmegaS=2*FS*tan(ws/2);
[k,Wn]=buttord(OmegaP,OmegaS,1,60,'s');
Wn=Wn/pi; %转换成数字滤波器的Wn
[b,a]=butter(k,Wn);
%freqs(b,a) %设计模拟的
%[bz,az]=bilinear(b,a,T); %映射为数字的
% 绘制结果
[H,f]=freqz(b,a,1024,'whole',FS);
%figure(8);
plot(f,abs(H)); grid;
title('Butterworth滤波器的频率响应');
xlabel('频率(Hz)');
ylabel('幅值');
figure(2);
[db, mag, pha, grd,w]=freqz_m(b,a);
plot(w*FS/2/pi,db); grid;
title('Butterworth滤波器的频率响应');
xlabel('频率(Hz)');
ylabel('幅值(dB)');
响应曲线如下
-----------------------------------
谢谢你们能来教我!
我的图和你前两个是一样的!但横坐标是0 100 200 300 400 500 600
你的0.1 那些是什么啊 ?这个是低通的吗?
我的图和你前两个是一样的!但横坐标是0 100 200 300 400 500 600
你的0.1 那些是什么啊 ?这个是低通的吗?
----------------------------------
你的图从0~1024是样点数，没有转成频率（包括你引用的帖子：http://www.chinavib.com/forum/thread-36508-1-1.html，它的图中给出的也是样点数），而我给的图，blackman滤波器中因你没有给出采样频率，用归正化频率表示；Butterworth滤波器中你给出了采样频率FS=1，横轴便是频率。ds2a图中0.1 是相对归一化频率（fs/2）的值；ds3a图中0.1 是0.1 Hz。
-------------------------------
Butterworth滤波器设计
FS=1;
%通带、阻带截止频率
Fl=0.08;Fh=0.15;
%频率预畸
wp=(Fl/FS)*2*pi; %临界频率采用角频率表示
ws=(Fh/FS)*2*pi; %临界频率采用角频率表示
OmegaP=2*FS*tan(wp/2);
OmegaS=2*FS*tan(ws/2);
[k,Wn]=buttord(OmegaP,OmegaS,1.1,60,'s');
FS是什么参数呀？？FL和FH 是如何确定的？？？
还有为什么要换成OmegaP,OmegaS这两个参数呢
我们现在在做这个课程设计，太恼火了！！哪位大哥知道，指导一下哈！！
FS=1;
%通带、阻带截止频率
Fl=0.08;Fh=0.15;
%频率预畸
wp=(Fl/FS)*2*pi; %临界频率采用角频率表示
ws=(Fh/FS)*2*pi; %临界频率采用角频率表示
OmegaP=2*FS*tan(wp/2);
OmegaS=2*FS*tan(ws/2);
[k,Wn]=buttord(OmegaP,OmegaS,1.1,60,'s');
FS是什么参数呀？？FL和FH 是如何确定的？？？
还有为什么要换成OmegaP,OmegaS这两个参数呢
我们现在在做这个课程设计，太恼火了！！哪位大哥知道，指导一下哈！！
--------------------------------
引用:
原帖由 tianma 于 2007-7-10 21:59 发表
Butterworth滤波器设计
FS=1;
%通带、阻带截止频率
Fl=0.08;Fh=0.15;
%频率预畸
wp=(Fl/FS)*2*pi; %临界频率采用角频率表示
ws=(Fh/FS)*2*pi; %临界频率采用角频率表示
OmegaP=2*FS*tan(wp/2);
OmegaS=2*FS*tan(ws/2);
[k,Wn]=buttord(OmegaP,OmegaS,1.1,60,'s');
FS是什么参数呀？？FL和FH 是如何确定的？？？
还有为什么要换成OmegaP,OmegaS这两个参数呢
我们现在在做这个课程设计，太恼火了！！哪位大哥知道，指导一下哈！！
这是用双线性Z变换的方法设计数字滤波器。FS是采样频率，FL和FH 是设计低通滤波器中的通带和阻带，根据具体要求来设定。因为用的方法是双线性Z变换，要对频率先进行预畸(wrap)，其原因可参看数字信号处理的书藉。
fs是采样率。fl，fh是根据你设计的滤波器的要求而定的，为什么要转换，你看一下滤波器设计的参考书，都会有的！！！
------------------------------
我问一个很弱智的问题，我设计出了滤波器但是不怎么会使用啊，谁能指教一下？
-----------------------------
引用:
原帖由 zhangnan3509 于 2007-8-28 23:11 发表
我问一个很弱智的问题，我设计出了滤波器但是不怎么会使用啊，谁能指教一下？
如果已设计出了滤波器，就是求出了数字滤波器的系数：a和b，那可以用filter函数对输入序列进行滤波了，设输入序列为x，输出序列为y，有
y=filter(b,a,x)
--------------------------------
y=filter(b,a,x)
--------------------------------
可能是我的仿真信号频率有问题，所以出的结果不正确 我再好好看看。
可我不知道为什么总得到这样的结果，滤波后的时频图是条直线啊
------------------------------
能把程序传上来，我看看。
---------------------------
回复 #15 songzy41 的帖子
clear;
f1 = 1776;
f2 = 1000;
fs = 4000;
n = 1:4000;
x = sin(2*pi*f1*n/fs);
x(1400:2900) = 0;
y = x + sin(2*pi*f2*n/fs);
Wp = 1500/2500; Ws = 1200/2500;
[n,Wn] = buttord(Wp,Ws,-0.05,3)
n=5
Wn=0.47985
[b,a]=butter(n,Wn)
[h,w]=freqz(b,a,100,4000)
figure(9)
subplot(2,1,1)
plot(w,abs(h));
%添加横向座标轴的标注
xlabel('频率(HZ)');
%添加纵向座标轴的标注
ylabel('幅值');
grid on
sf=filter(b,a,y); %叠加函数x经过低通滤波器以后的新函数
subplot(212);
plot(n/fs,sf); %绘制叠加函数x经过低通滤波器以后的时域图形
xlabel('时间 (seconds)');
ylabel('幅度');
不知道为什么出来的是低通，而且滤波后还是直线？
f1 = 1776;
f2 = 1000;
fs = 4000;
n = 1:4000;
x = sin(2*pi*f1*n/fs);
x(1400:2900) = 0;
y = x + sin(2*pi*f2*n/fs);
Wp = 1500/2500; Ws = 1200/2500;
[n,Wn] = buttord(Wp,Ws,-0.05,3)
n=5
Wn=0.47985
[b,a]=butter(n,Wn)
[h,w]=freqz(b,a,100,4000)
figure(9)
subplot(2,1,1)
plot(w,abs(h));
%添加横向座标轴的标注
xlabel('频率(HZ)');
%添加纵向座标轴的标注
ylabel('幅值');
grid on
sf=filter(b,a,y); %叠加函数x经过低通滤波器以后的新函数
subplot(212);
plot(n/fs,sf); %绘制叠加函数x经过低通滤波器以后的时域图形
xlabel('时间 (seconds)');
ylabel('幅度');
不知道为什么出来的是低通，而且滤波后还是直线？
-------------------------------
主要的错误是在滤波器的设计中。在 buttord(Wp,Ws,Rp,Rs)中，Wp是通带的频率、Ws是阻带的频率、Rp是通带内波纹系数需小于的值、Rs是阻带内衰减系数需大于的值。Rp和Rs都表示的分贝值，都为正值。若设计低通滤波器，有Wp<Ws。同时Wp和Ws都是归一化频率，设滤波器的实际通带频率为fp和实际阻带频率为fc，则Wp=fp/(fs/2)、Ws=fc/(fs/2)。
从版主的图来看是设计低通滤波器，但参数中Wp = 1500/2500; Ws = 1200/2500; Wp>Ws(?)设计高通可Wp>Ws。
又Wp=fp/(fs/2)，应该除2000，怎么会除2500？
所以设计的滤波器完全不能满足要求，滤波器输出也就怪怪的。我把程序改成：
clear;
f1 = 1776;
f2 = 1000;
fs = 4000;
n = 1:4000;
x = sin(2*pi*f1*n/fs);
x(1400:2900) = 0;
y = x + sin(2*pi*f2*n/fs);
Wp = 1200/2000; Ws = 1500/2000;
[n,Wn] = buttord(Wp,Ws,0.05,3)
%n=5;
%Wn=0.47985;
[b,a]=butter(n,Wn);
[h,w]=freqz(b,a);
figure(9)
subplot(2,1,1)
plot(fs*w/2/pi,abs(h));
%添加横向座标轴的标注
xlabel('频率(HZ)');
%添加纵向座标轴的标注
ylabel('幅值');
grid on
%break
sf=filter(b,a,y); %叠加函数x经过低通滤波器以后的新函数
subplot(212);
plot(sf); %绘制叠加函数x经过低通滤波器以后的时域图形
xlabel('时间 (seconds)');
ylabel('幅度');
作出的图如下。
从版主的图来看是设计低通滤波器，但参数中Wp = 1500/2500; Ws = 1200/2500; Wp>Ws(?)设计高通可Wp>Ws。
又Wp=fp/(fs/2)，应该除2000，怎么会除2500？
所以设计的滤波器完全不能满足要求，滤波器输出也就怪怪的。我把程序改成：
clear;
f1 = 1776;
f2 = 1000;
fs = 4000;
n = 1:4000;
x = sin(2*pi*f1*n/fs);
x(1400:2900) = 0;
y = x + sin(2*pi*f2*n/fs);
Wp = 1200/2000; Ws = 1500/2000;
[n,Wn] = buttord(Wp,Ws,0.05,3)
%n=5;
%Wn=0.47985;
[b,a]=butter(n,Wn);
[h,w]=freqz(b,a);
figure(9)
subplot(2,1,1)
plot(fs*w/2/pi,abs(h));
%添加横向座标轴的标注
xlabel('频率(HZ)');
%添加纵向座标轴的标注
ylabel('幅值');
grid on
%break
sf=filter(b,a,y); %叠加函数x经过低通滤波器以后的新函数
subplot(212);
plot(sf); %绘制叠加函数x经过低通滤波器以后的时域图形
xlabel('时间 (seconds)');
ylabel('幅度');
作出的图如下。
附件
--------------------------------
采样频率原来不是4000，而是5000，修改了一次忘了改回来，
。
另外一句是plot(fs*w/2/pi,abs(h))修改了这句plot(w,abs(h))，w返回是弧度，因此要做转化，我没看仔细help
。另外一句是plot(fs*w/2/pi,abs(h))修改了这句plot(w,abs(h))，w返回是弧度，因此要做转化，我没看仔细help
为什么改变了wp,ws的值一样还是低通的？比如把低通1200/2000，1500/2000这两个值对调一下，应该就是高通了吧
------------------------------------
再把[b,a]=butter(n,Wn);
改成b,a]=butter(n,Wn,'high');就可以了
其实这样编程来设计滤波器太麻烦了，可以用FDATool来设计。
它完全是图形化的，能设计各类滤波器，并且可以很方便地导出滤波器系数。
改成b,a]=butter(n,Wn,'high');就可以了
其实这样编程来设计滤波器太麻烦了，可以用FDATool来设计。
它完全是图形化的，能设计各类滤波器，并且可以很方便地导出滤波器系数。
-----------------------------------
……
-------------------------------
freqz_m是什么函数？是不是freqz？你写错函数名了吧
-------------------------------
有freqz_m这个函数，它对freqz作了一些改进。该函数出自“数字信号处理--使用MATLAB” （维纳。恩格尔编著），该函数为
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);
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);
-----------------------------------
- 上一篇： 低通滤波的效果不理想
- 下一篇： fft采用不同分析点数为什么幅值会不一样呢？ 频谱校正
EE直播间
更多
-
提升毫米波信号测试精度 直播时间： 12月18日 14:00 -
EE Talk主题专访系列直播-对话：释放 Wi-Fi 7 在高带宽应用中的技术潜力 直播时间： 12月19日 10:00
关闭
站长推荐
/2 
/2
文章评论（0条评论）
登录后参与讨论