请教关于低通滤波器设计的问题

zhangnan3509

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('幅度');
不知道为什么出来的是低通，而且滤波后还是直线？

songzy41

回复 #16 zhangnan3509

主要的错误是在滤波器的设计中。在 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('幅度');

作出的图如下。

zhangnan3509

采样频率原来不是4000，而是5000，修改了一次忘了改回来，:loveliness: 。
另外一句是plot(fs*w/2/pi,abs(h))修改了这句plot(w,abs(h))，w返回是弧度，因此要做转化，我没看仔细help

zhangnan3509

为什么改变了wp,ws的值一样还是低通的？比如把低通1200/2000，1500/2000这两个值对调一下，应该就是高通了吧

破凰

再把[b,a]=butter(n,Wn);
改成b,a]=butter(n,Wn,'high');就可以了
其实这样编程来设计滤波器太麻烦了，可以用FDATool来设计。
它完全是图形化的，能设计各类滤波器，并且可以很方便地导出滤波器系数。

弯弓射大雕

老师好，我是土木的，对信号处理一点不懂，现在想对采集的数据进行低通滤波，不知用什么滤波程序呢，怎么用啊？我采集的数据点数为9600，采样间隔为0.05s，想把0.02Hz以上的频率滤掉。谢谢了

[ 本帖最后由 zhlong 于 2007-9-15 12:38 编辑 ]

zpact

我有一个很弱智的问题，大家别打我。我在运行设计好的butterworth滤波器。 系统显示??? Undefined command/function 'buttord'.
我用help buttord; help butter; help chevby1;  都说这个function的文件都找不到。我知道这是个很初级的问题。但希望大家解答下。

[ 本帖最后由 zpact 于 2007-10-30 08:34 编辑 ]

zhangnan3509

你的MATLAB里面没有这个函数，是不是版本的问题啊？

zpact

我是7.1的版本，没理由，没有啊，我在学校7.4的版本上也试验了，没有。Cheby1,Ellip,Butter,Butterord.

zhangnan3509

which buttord
D:\Program Files\MATLAB71\toolbox\signal\signal\buttord.m
which  butter
D:\Program Files\MATLAB71\toolbox\signal\signal\butter.m

which Cheby1
D:\Program Files\MATLAB71\toolbox\signal\signal\cheby1.m
which Ellip
D:\Program Files\MATLAB71\toolbox\signal\signal\ellip.m

我的也是7.1，为什么我都有呢？是不是文件损坏了？

zpact

众所周知的原因，发言屏蔽了
                            --------zhangnan3509

[ 本帖最后由 zhangnan3509 于 2007-10-30 10:14 编辑 ]

zhangnan3509

以后这些信息请在本版置顶贴求助区里面发帖求助，为了防止有水贴出现在讨论区，谢谢合作

zpact

虽然被屏蔽了，但还是很开心，谢谢张楠同学，因为我在国外，所以这个问题实在不好解决，多谢你解我燃眉之急。

zpact

:'( 可能是matlab公司有什么防copy技术吧，我复制了，但是还是没有办法使用这些function。我还是想办法下载个好版本吧。

zpact

>> run('C:\Program Files\MATLAB71\toolbox\signal\signal\buttord.m')
??? Error using ==> run
Error using ==> evalin
Undefined function or variable 'buttord'.