butter函数：巴特沃斯模拟和数字滤波器设计

详细解释如下：

butter Butterworth digital and analog filter design.
   [B,A] = butter(N,Wn) designs an Nth order lowpass digital
   Butterworth filter and returns the filter coefficients in length
   N+1 vectors B (numerator) and A (denominator). The coefficients
   are listed in descending powers of z. The cutoff frequency
   Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to
   half the sample rate.

   If Wn is a two-element vector, Wn = [W1 W2], butter returns an
   order 2N bandpass filter with passband  W1 < W < W2.
   [B,A] = butter(N,Wn,'high') designs a highpass filter.
   [B,A] = butter(N,Wn,'low') designs a lowpass filter.
   [B,A] = butter(N,Wn,'stop') is a bandstop filter if Wn = [W1 W2].
   When used with three left-hand arguments, as in
   [Z,P,K] = butter(...), the zeros and poles are returned in
   length N column vectors Z and P, and the gain in scalar K.

   When used with four left-hand arguments, as in
   [A,B,C,D] = butter(...), state-space matrices are returned.

   butter(N,Wn,'s'), butter(N,Wn,'high','s') and butter(N,Wn,'stop','s')
   design analog Butterworth filters.  In this case, Wn is in [rad/s]
   and it can be greater than 1.0

需要注意的是，Wn不能大于1，Wn = 1相当于fs/2。

例1：

t = 0:0.001:1;
x = sin(t*2*pi*200)+ sin(t*2*pi);

[b,a] = butter(10,0.2);
figure(1);
freqz(b,a);

y = filter(b,a,x);
figure(2)
subplot(2,1,1);
plot(x);

subplot(2,1,2);
plot(y);

结果：

http://s9/middle/84024a4a4cb3618adf7e8&amp;690

http://s5/middle/84024a4a4cb3618c36ed4&amp;690

例2：生成一个采样频率为2000，由50Hz、120Hz、200Hz、正弦信号及噪声信号组成的数字信号，滤掉小于150Hz的频率，计算并显示滤波前后的原始数据波形以及功率谱密度

下面是代码，其中[Y,f] = Spectrum_Calc(y,Fs);函数的具体代码见http://blog.sina.com.cn/s/blog_84024a4a01015rez.html

close all
clear all
clc
Fs=2000;
Fa=Fs/2;
N=512;
t=0:1/Fs:1;
x=2*sin(t*2*pi*50)+sin(t*2*pi*120)+3*sin(t*2*pi*200)+randn(1,length(t));

% [b,a]=butter(10,[100/Fa 150/Fa],'stop');
[b,a]=butter(10,150/Fa,'high');
figure(1);freqz(b,a);
y=filter(b,a,x);

figure(2);
subplot(2,1,1);plot(x);
subplot(2,1,2);plot(y);

figure(3);
subplot(2,1,1);
[Y,f] = Spectrum_Calc(x,Fs);
P=Y.*conj(Y)/N;  
plot(f,Y);

figure(3);
subplot(2,1,2);
[Y,f] = Spectrum_Calc(y,Fs);
P=Y.*conj(Y)/N;  
plot(f,Y);

运行结果：

http://s7/middle/84024a4a07ab8b4d35f06&amp;690

http://s12/middle/84024a4a4cb371047d00b&amp;690

fourier变换的幅频图

http://s1/middle/84024a4a4cb3710563af0&amp;690

功率谱密度图，这里需要将上面的plot(f,Y);换为plot(f,P);

http://s4/middle/84024a4a4cb37106009f3&amp;690

从功率谱和幅频图上看，效果是很好的。

另外我还尝试了带阻滤波器，但是效果不好，无法得到预期图像。