关于用MATLAB设计对信号进行频谱分析和滤波处理的程序

完整的程序

%写上标题

%设计低通滤波器：

[N,Wc]=buttord()

%估算得到Butterworth低通滤波器的最小阶数N和3dB截止频率Wc

[a,b]=butter(N,Wc) %设计Butterworth低通滤波器

[h,f]=freqz() %求数字低通滤波器的频率响应

figure(2) % 打开窗口2

subplot(221) %图形显示分割窗口

plot(f,abs(h)) %绘制Butterworth低通滤波器的幅频响应图

title(巴氏低通滤波器'')

grid %绘制带网格的图像

sf=filter(a,b,s) %叠加函数S经过低通滤波器以后的新函数

subplot(222)

plot(t,sf) %绘制叠加函数S经过低通滤波器以后的时域图形

xlabel('时间 (seconds)')

ylabel('时间按幅度')

SF=fft(sf,256) %对叠加函数S经过低通滤波器以后的新函数进行256点的基—2快速傅立叶变换

w= %新信号角频率

subplot(223)

plot())%绘制叠加函数S经过低通滤波器以后的频谱图

title('低通滤波后的频谱图')

%设计高通滤波器

[N,Wc]=buttord()

%估算得到Butterworth高通滤波器的最小阶数N和3dB截止频率Wc

[a,b]=butter(N,Wc,'high')%设计Butterworth高通滤波器

[h,f]=freqz()%求数字高通滤波器的频率响应

figure(3)

subplot(221)

plot())%绘制Butterworth高通滤波器的幅频响应图

title('巴氏高通滤波器')

grid%绘制带网格的图像

sf=filter()%叠加函数S经过高通滤波器以后的新函数

subplot(222)

plot(t,sf)%绘制叠加函数S经过高通滤波器以后的时域图形

xlabel('Time(seconds)')

ylabel('Time waveform')

w%新信号角频率

subplot(223)

plot())%绘制叠加函数S经过高通滤波器以后的频谱图

title('高通滤波后的频谱图')

%设计带通滤波器

[N,Wc]=buttord([)

%估算得到Butterworth带通滤波器的最小阶数N和3dB截止频率Wc

[a,b]=butter(N,Wc)%设计Butterworth带通滤波器

[h,f]=freqz()%求数字带通滤波器的频率响应

figure(4)

subplot(221)

plot(f,abs(h))%绘制Butterworth带通滤波器的幅频响应图

title('butter bandpass filter')

grid%绘制带网格的图像

sf=filter(a,b,s)%叠加函数S经过带通滤波器以后的新函数

subplot(222)

plot(t,sf)%绘制叠加函数S经过带通滤波器以后的时域图形

xlabel('Time(seconds)')

ylabel('Time waveform')

SF=fft()%对叠加函数S经过带通滤波器以后的新函数进行256点的基—2快速傅立叶变换

w=( %新信号角频率

subplot(223)

plot('))%绘制叠加函数S经过带通滤波器以后的频谱图

title('带通滤波后的频谱图')

采样没什么，就是产生一个连续的（实际还是数字信号），实际上就是再进行一下抽取。

变化就用fft函数。

滤波器设计有专门的函数来实现，IIR的有巴特沃斯、切比雪夫、椭圆等等。FIR可以直接在频域设计，应该也有专门的函数，忘了。高通就是用1减去低通，带通就是高通加低通减1，当然这是比较投机的方法，数字信号处理在FIR设计里有专门讲几种滤波器的设计。

用buttord和buffer得到了拉普拉斯变换的分子分母多项式系数a，b,

假设信号是x，则就用y=filter(b,a,x)

例如：

设计一个高通滤波器，并检验它的性能

采样率为10kHZ

阻带边缘为1.5Khz,衰减为40bB

通带边缘为2kHz，波纹为3Db

>>Fs=1e4

>>fs=1.5e3

>>fp=2e3

>>As=40

>>Rp=3

>>wp=2*fp/Fs

>>ws=2*fs/Fs

>>[N,wn]=cheb2ord(wp,ws,Rp,As)

>>[b,a]=cheby2(N,As,wn,'high')

>>[db,mag,pha,grd,w]=freqz_m(b,a)

>>subplot(2,2,1)plot(w/pi,mag)

>>axis([0,1,0,1])

>>setX([0 0.3 0.4 1])

>>setY([0.01 0.7279 1])

>>title('Magnitude Response')

>>subplot(2,2,2)plot(w/pi,db)

>>axis([0 1 -70 0])

>>setX([0 0.3 0.4 1])

>>setY([-40 -2.7589])

>>title('Magnitude Response in dB')

然后给你一个信号x=cos(0.2*pi*n)

>>n=0:200

>>x=cos(0.6*pi*n)

>>y=filter(b,a,x)

>>subplot(2,2,3)plot(n,x)

>>subplot(2,2,4)plot(n,y)

>>x1=fft(x,201)

>>x11=abs(x1)

>>subplot(2,2,1)stem(n,x11)

>>y1=fft(y,201)

>>y11=abs(y1)

>>subplot(2,2,2)stem(n,y11)

>>setX([0 60 140 201])

>>title('FFT of y')

>>subplot(2,2,1)stem(n,x11)

>>setX([0 60 140 201])

>>title('FFT of x')

>>g=x11-y11

>>subplot(2,2,3)stem(n,g)

你运行一下，就可以看到结果，刚好把此信号滤掉

---