DSP实验报告-南邮.doc

上机实验内容:
数组的加、减、乘、除和乘方运算。输入A=[1 2 3 4]，B=[3 4 5 6]，求C=A+B;
D=A-B;E=A.*B;F=A./B;G=A.^B;并用Stem语句画出A、B、C、D、E、F、G。
(A)clear;
A=[1 2 3 4];
stem(A);
A:

(B)
B=[3,4,5,6];
stem(B);

B:

(C)
A=[1 2 3 4];
B=[3,4,5,6];
C=A+B

C =

4 6 8 10
A=[1 2 3 4];
B=[3,4,5,6];
C=A+B

C(1)

ans =

4

>>
>> A=[1 2 3 4];
B=[3,4,5,6];
C=A+B;
C(1:3)

ans =

4 6 8

>>
>> A=[1 2 3 4];
%stem(A);
B=[3,4,5,6];
%stem(B);
C=A+B;
%C(1);
%C(1:3);
C(0)
??? Subscript indices must either be real positive integers or logicals.

>>
A=[1 2 3 4];
B=[3,4,5,6];
C=A+B;
stem(C);

A=[1 2 3 4];
B=[3,4,5,6];
C=A+B;
n=0:1:3;
n;
stem(n,C);

C =

4 6 8 10

>>
(D)

A=[1 2 3 4];
B=[3,4,5,6];
D=A-B;
stem(D);

D:

D =

-2 -2 -2 -2

>>
(E)
A=[1 2 3 4];
B=[3,4,5,6];
C=A+B;
E=A.*B;
stem(E);

E:

>>E=A.*B
E =

3 8 15 24
(F)

A=[1 2 3 4];
B=[3,4,5,6];
F=A./B;
stem(F);

F:

>> F=A./B

F =

0.3333 0.5000 0.6000 0.6667
(G)
A=[1 2 3 4];
B=[3,4,5,6];
G=A.^B;
stem(G);

>> G=A.^B

G =

1 16 243 4096
(2) 用MATLAB实现下列序列:

①
clear;
n=0:1:15;
x1=0.8.^n;
n=0:1:15;
n;
stem(n,x1);

②
n=0:1:15;
a=(0.2+3*j)*n;
x2=exp(a);
n=0:1:15;
n;
figure;
stem(n,x2);

③
n=0:1:15;
x3=3*cos(0.125*pi*n+0.2*pi)+2*sin(0.25*pi*n+0.1*pi);
figure;
stem(n,x3);

d) 将(c)中的x(n)扩展为以16为周期的函数 , 绘出四个周期。
④
n=0:1:63;
x4=3*cos(0.125*pi*rem(n,16)+0.2*pi)+2*sin(0.25*pi*rem(n,16)+0.1*pi);
figure;
stem(n,x4);

e) 将(c)中的x(n)扩展为以10为周期的函数，绘出四个周期。
⑤
n=0:1:39;
x5=3*cos(0.125*pi*rem(n,10)+0.2*pi)+2*sin(0.25*pi*rem(n,10)+0.1*pi);
figure;
stem(n,x5);

(3) 产生并绘出下列序列的样本:

n=0:3;
x=[1 -1 3 5 ]
x1=circshift(x,[0,-2]);
x2=circshift(x,[0,1]);
x3=2*x1-x2-2*x;
stem(n,x3);

x=[1 -1 3 5];
x1=zeros(1,4);
x2=zeros(1,4);
for k=1:5
for n=0:3
x3=circshift(x,[0,-k])
t=n.*x3;
x1=x1+t;
end
x2=x2+x1;
end
stem(x2);

(4) 绘出下列时间函数的图形，对x轴、y轴以及图形上方均须加上适当的标注

t=0:0.001:10;
x=sin(2*pi*t);
plot(t,x)
title('x=sin(2*pi*t)');
xlabel('x');
ylabel('t');

t=0:0.01:10;
x=sin(2*pi*t);
plot(t,x);
plot(t,x,'r_');
title('sin');
xlabel('t');
ylabel('x(t)');
title('sin');

t=0:0.01:4;
x1=cos(100*pi*t);
x2=sin(pi*t);
x=x1.*x2;
plot(t,x)
title('x=cos(100*pi*t)*sin(pi*t)');
xlabel('x');
ylabel('t');

(5) 编写函数 实现，绘出该函数的图形，起点为n1，终点为n2。
function stepshift(n0,n1,n2); %单位阶跃序列,n0为时移量
n=n1:n0-1; %n1、n2为序列的起止序列号
nn=length(n);
x=zeros(1,nn); %n0前信号赋值为0
stem(n,x,'fill') %绘出n1~n0-1的波形(0值)
hold on
k=n0:n2;
kk=length(k);
x=ones(1,kk); %n0后信号赋值为1
stem(k,x,'fill') %绘出n1~n0-1的波形(1值)
hold off
axis([n1,n2,0,1.1])
title('单位阶跃序列')
运行程序后，在Commond Windows中输入stepshift(6,-3,24)
则得到平移后序列图形如下：

(6) 给定一因果系统 求出并绘制H(z)的幅频响应与相频响应。
clear all;
k=64;
b=[1 sqrt(2)];
a=[1 -0.67 0.9];
w=0:pi/k:pi;
h=freqz(b,a,w);
subplot(2,1,1);
plot(abs(h))
grid on
subplot(2,1,2);
plot(angle(h))
grid on

(7) 计算序列和序列的离散卷积，并作图表示卷积结果。
A=[8 -2 -1 2 3 ]
B=[2 3 -1 -1];
C=conv(A ,B)
plot(C)

(8) 求以下差分方程所描述系统的单位脉冲响应

clear all;
N=50;
a=[1 -2];
b=[1 0.1 -0.06];
x1=[1 zeros(1,N-1)];
n=0:1:N-1;
h=filter(a,b,x1);
stem(n,h)
axis([-1 53 -2.5 1.2])

三：思考题

（1）对于有限长序列，如何用MATLAB计算其DTFT?

fs=1000;
t=0:1/fs:0.6;
f1=100;
f2=300;
x=sin(2*pi*f1*t)+sin(2*pi*f2*t);
subplot(711)
plot(x);
title(‘f1(100hz)\f2(300hz)的正弦信号，初相0’)
xlabel(‘序列（n）’)
grid on

number=512
y=fft(x,number);
n=0:length(y)-1;
f=fs*n/length(y);
subplot(713)
plot(f,abs(y));
title('f1\f2的正弦信号的fft（512点）')
xlabel('频率hz')
grid on

x=x+randn(1,length(x));
subplot(715)
plot(x);
title('原f1\f2的正弦信号（含随机噪声）')
xlabel('序列（n）')
grid on

y=fft(x,number);
n=0:length(y)-1;
f=fs*n/length(y);
subplot(717)
plot(f,abs(y));
title('原f1\f2的正弦信号（含随机噪声）的fft（512点）')
xlabel('频率hz')
grid on

(3) 对于由两个子系统级联或并联的系统，如何用MATLAB计算它们的幅频响应与相频响应?
答：级联转换为直接型：cas2dir 并联转换为直接型：par2dir 然后用freqz（）就行了。
实验二 快速傅里叶变换(FFT)及其应用
一、实验目的
(1) 在理论学习的基础上，通过本实验，加深对FFT的理解，熟悉MATLAB中的有关函数。
(2) 应用FFT对典型信号进行频谱分析。
(3) 了解应用FFT进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT。
(4) 应用FFT实现序列的线性卷积和相关。
二、实验内容
实验中用到的信号序列:
a) 高斯序列

b) 衰减正弦序列

c) 三角波序列

d) 反三角波序列

上机实验内容:

方法二：

n=0:15;
p=8;
q=2;
xa=exp(-(n-p).^2/q);
subplot(3,2,1);
stem(n,xa,'.');
title('q=2高斯序列时域特性');
subplot(3,2,2);
stem(fft(xa),'.');
title('q=2高斯序列幅频特性');

n=0:15;
p=8;
q=4;
xa=exp(-(n-p).^2/q);
subplot(3,2,3);
stem(n,xa,'.');
title('q=4高斯序列时域特性');
subplot(3,2,4);
stem(fft(xa),'.');
title('q=4高斯序列幅频特性');

n=0:15;
p=8;
q=8;
xa=exp(-(n-p).^2/q);
subplot(3,2,5);
stem(n,xa,'.');
title('q=8高斯序列时域特性');
subplot(3,2,6);
stem(fft(xa),'.');
title('q=8高斯序列幅频特性');

n=0:15;
p=8;
q=8;
xa=exp(-(n-p).^2/q);
subplot(3,2,1);
stem(n,xa,'.');
title('p=8高斯序列时域特性');
subplot(3,2,2);
stem(fft(xa),'.');
title('p=8高斯序列幅频特性');

n=0:15;
p=13;
q=8;
xa=exp(-(n-p).^2/q);
subplot(3,2,3);
stem(n,xa,'.');
title('p=13高斯序列时域特性');
subplot(3,2,4);
stem(fft(xa),'.');
title('p=13高斯序列幅频特性');

n=0:15;
p=14;
q=8;
xa=exp(-(n-p).^2/q);
subplot(3,2,5);
stem(n,xa,'.');
title('p=14高斯序列时域特性');
subplot(3,2,6);
stem(fft(xa),'.');
title('p=14高斯序列幅频特性');

n=0:15;
p=8;
q=8;
xa=exp(-(n-p).^2/q);
subplot(3,2,1);
stem(n,xa,'.');
title('p=8高斯序列时域特性');
subplot(3,2,2);
plot(abs(fft(xa)));
title(' p=8高斯序列幅频特性');

n=0:15;
p=13;
q=8;
xa=exp(-(n-p).^2/q);
subplot(3,2,3);
stem(n,xa,'.');
title(' p=13高斯序列时域特性');
subplot(3,2,4);
plot(abs(fft(xa)));
title(' p=13高斯序列幅频特性');

n=0:15;
p=14;
q=8;
xa=exp(-(n-p).^2/q);
subplot(3,2,5);
stem(n,xa,'.');
title(' p=14高斯序列时域特性');
subplot(3,2,6);
plot(abs(fft(xa)));
title(' p=14高斯序列幅频特性');

(2)
n=0:1:15;
A=1;
f=0.0625;
a=0.1;
x=exp(-a*n).*sin(2*pi*f*n);
subplot(3,2,1);
stem(n,x,'.');
axis([0 15 -2 2]);
title('f=0.0625衰减正弦序列时域特性');
subplot(3,2,2);
stem(fft(x),'.');
axis([0 15 -5 2]);
title('f=0.0625衰减正弦序列幅频特性');

n=0:1:15;
A=1;
f=0.4375;
a=0.1;
x=exp(-a*n).*sin(2*pi*f*n);
subplot(3,2,3);
stem(n,x,'.');
axis([0 15 -2 2]);
title('f=0.4375衰减正弦序列时域特性');
subplot(3,2,4);
stem(fft(x),'.');
axis([0 15 -5 2]);
title('f=0.4375衰减正弦序列幅频特性');

n=0:1:15;
A=1;
f=0.5625;
a=0.1;
x=exp(-a*n).*sin(2*pi*f*n);
subplot(3,2,5);
stem(n,x,'.');
axis([0 15 -2 2]);
title('f=0.5625衰减正弦序列时域特性');
subplot(3,2,6);
stem(fft(x),'.');
axis([0 15 -5 2]);
title('f=0.5625衰减正弦序列幅频特性');

n=0:1:15;
A=1;
f=0.0625;
a=0.1;
x=exp(-a*n).*sin(2*pi*f*n);
subplot(3,2,1);
stem(n,x,'.');
axis([0 15 -2 2]);
title('f=0.0625衰减正弦序列时域特性');
subplot(3,2,2);
plot(abs(fft(x)));
axis([0 20 0 5]);
title('f=0.0625衰减正弦序列幅频特性');

n=0:1:15;
A=1;
f=0.4375;
a=0.1;
x=exp(-a*n).*sin(2*pi*f*n);
subplot(3,2,3);
stem(n,x,'.');
axis([0 15 -2 2]);
title('f=0.4375衰减正弦序列时域特性');
subplot(3,2,4);
plot(abs(fft(x)));
axis([0 20 0 5]);
title('f=0.4375衰减正弦序列幅频特性');

n=0:1:15;
A=1;
f=0.5625;
a=0.1;
x=exp(-a*n).*sin(2*pi*f*n);
subplot(3,2,5);
stem(n,x,'.');
axis([0 15 -2 2]);
title('f=0.5625衰减正弦序列时域特性');
subplot(3,2,6);
plot(abs(fft(x)));
axis([0 20 0 5]);
title('f=0.5625衰减正弦序列幅频特性');

(3)

clear;
for n=1:4;
xc(n)=n-1;
end
for n=5:8;
xc(n)=8-(n-1);
end
m=0:7;
subplot(2,3,1);
stem(m,xc,'.');
title('三角波序列');

Xc=fft(xc,8);
k=0:1:7;
subplot(2,3,2);
stem(k,abs(Xc),'.');
title('三角波序列8点FFT');
Xc=fft(xc,32);
k=0:1:31;
subplot(2,3,3);
stem(k,abs(Xc),'.');
title('三角波系列32点FFT');

clear;
for n=1:4;
xc(n)=4-(n-1);
end
for n=5:8;
xc(n)=(n-1)-4;
end
m=0:7;
subplot(2,3,4);
stem(m,xc,'.');
title('反三角波序列');

Xc=fft(xc,8);
k=0:1:7;
subplot(2,3,5);
stem(k,abs(Xc),'.');
title('反三角波序列8点FFT');
Xc=fft(xc,32);
k=0:1:31;
subplot(2,3,6);
stem(k,abs(Xc),'.');
title('反三角波序列32点FFT');

(4) N=16,△f=1/16时，其频谱：

N=16;
n=0:1:N-1;
x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+1/16)*n);
subplot(2,1,1);
plot(n,x,'.');
axis([0 15 -2 2]);
grid;
title('连续时间信号');
subplot(2,1,2);
plot(abs(fft(x)));
axis([0 20 0 10]);
grid;
title('连续时间信号频谱');

N=16,△f=1/64时，其频谱：

N=16;
n=0:1:N-1;
x=sin(2*pi*0.125*n)+cos(2*pi*(0.125+1/64)*n);
subplot(2,1,1);
plot(n,x,'.');
axis([0 15 -2 2]);
grid;
title('连续时间信号');
subplot(2,1,2);
plot(abs(fft(x)));
axis([0 20 0 10]);
grid;
title('连续时间信号频谱');

N=128,△f=1/16时，其频谱：

N=128,△f=1/64时，其频谱：

(5)循环卷积：
n=0:15;
p=8;
q=2;
xa=exp(-(n-p).^2/q);
subplot(2,3,1);
stem(n,xa,'.');
title('xa波形');
Xa=fft(xa,16);
subplot(2,3,4);
stem(abs(Xa),'.');
title('Xa(k)=FFT[xa(n)]的波形 ');

n=0:1:15;
A=1;
f=0.0625;
a=0.1;
xb=exp(-a*n).*sin(2*pi*f*n);
subplot(2,3,2);
stem(n,xb,'.');
title('xb波形');
Xb=fft(xb,16);
subplot(2,3,5);
stem(abs(Xb),'.');
title('Xb(k)=FFT[xb(n)]的波形 ');

Y=Xa.*Xb;
y=ifft(Y,16);
subplot(2,3,3);
stem(abs(Y),'.');
title('Y(k)=Xa(k)Xb(k)的波形);

subplot(2,3,6);
stem(y,'.');
title('y=IFFT[Y(k)]的波形 ');

线性卷积：
n=0:15;
p=8;
q=2;
xa=exp(-(n-p).^2/q);
subplot(2,3,1);
stem(n,xa,'.');
title('xa波形');
Xa=fft(xa,16);
subplot(2,3,4);
stem(abs(Xa),'.');
title('Xa(k)=FFT[xa(n)]的波形 ');

n=0:1:15;
f=0.0625;
a=0.1;
xb=exp(-a*n).*sin(2*pi*f*n);
subplot(2,3,2);
stem(n,xb,'.');
title('xb波形');
Xb=fft(xb,16);
subplot(2,3,5);
stem(abs(Xb),'.');
title('Xb(k)=FFT[xb(n)]的波形 ');
N1=length(xa);
N2=length(xb);
N=N1+N2-1;
xa=[xa zeros(1,N2-1)];
xb=[xb zeros(1,N1-1)];
n=1:N;
k=n;
Xa=fft(xa);
Xb=fft(xb);
Y=Xa.*Xb;
subplot(2,3,3);
stem(abs(Y),'.');
title('Y(k)=Xa(k).*Xb(k)的波形');
y=ifft(Y(k));
subplot(2,3,6);
stem(y,'.');
title('y=ifft[Y(k)]的波形 ');

六．
clear;

xe=rand(1,512);
for i=1:4
n(i)=i-1;
xc(i)=n(i);
end
for i=5:8
n(i)=i-1;
xc(i)=8-n(i);
end

%重叠相加法
yn=zeros(1,519);
for j=0:7
xj=xe(64*j+1:64*(j+1));

xak=fft(xj,71);
xck=fft(xc,71);
yn1=ifft(xak.*xck);

temp=zeros(1,519);

temp(64*j+1:64*j+71)=yn1;

yn=yn+temp;
end;
n=0:518;
figure(1)
subplot(211);
stem(n,yn);
xlabel('n');ylabel('y(n)');
title('重叠相加法:xc(n)与xe(n)的线性卷积的时域波形');
subplot(212);
stem(n,abs(fft(yn)));
xlabel('k');ylabel('Y(k)');
axis([0,600,0,300]);
title('重叠相加法:xc(n)与xe(n)的线性卷积的幅频特性');

%重叠保留法
k=1:7;
xe1=k-k;
xe_1=[xe1,xe];
yn_1=zeros(1,519);
for j=0:7
xj_1=xe_1(64*j+1:64*j+71);
xak_1=fft(xj_1);
xck_1=fft(xc,71);
yn1_1=ifft(xak_1.*xck_1);
temp=zeros(1,519);
temp(64*j+1:64*j+64)=yn1_1(8:71);
yn_1=yn_1+temp;

end;

n=0:518;
figure(2)
subplot(211);
stem(n,yn_1);
xlabel('n');ylabel('y(n)');
title(' 重叠保留法:xc(n)与xe(n)的线性卷积的时域波形');
subplot(212);
stem(n,abs(fft(yn_1)));
xlabel('k');ylabel('Y(k)');
axis([0,600,0,300]);
title('重叠保留法:xc(n)与xe(n)的线性卷积的幅频特性');

七．
clear;

N=16;
L=16;
n=0:N-1;
for i=1:N
n(i)=i-1;
xa(i)=exp(-(n(i)-8).*(n(i)-8)/2);
xb(i)=exp(-0.1*n(i))*sin(2*pi*0.0625*n(i));
end
Xa=fft(xa,2*N);
Xb=fft(xb,2*N);

rm=ifft(conj(Xa).*Xb);
rm=[rm(N+2:2*N) rm(1:N)];
m=[(-N+1):(N-1)];
subplot(221)
stem(m,rm);
title('线性相关rxy');

rm=ifft(conj(Xb).*Xa);
rm=[rm(N+2:2*N) rm(1:N)];
m=[(-N+1):(N-1)];
subplot(222)
stem(m,rm);
title('线性相关ryx');

Xa16=fft(xa,N);
Xb16=fft(xb,N);

rm=real(ifft(conj(Xa16).*Xb16));
subplot(223)
stem(rm);
title('循环相关rxy');

rm=real(ifft(conj(Xb16).*Xa16));
subplot(224)
stem(rm);
title('循环相关ryx');

八．
clear;

N=16;
L=16;
n=0:N-1;
for i=1:N
n(i)=i-1;
xa(i)=exp(-(n(i)-8).*(n(i)-8)/2);
xb(i)=exp(-0.1*n(i))*sin(2*pi*0.0625*n(i));
end
Xa=fft(xa,2*N);
Xb=fft(xb,2*N);

rm=real(ifft(conj(Xa).*Xa));
rm=[rm(N+2:2*N) rm(1:N)];
m=[(-N+1):(N-1)];
subplot(211)
stem(m,rm);
title('xa自相关');

rm=ifft(conj(Xb).*Xb);
rm=[rm(N+2:2*N) rm(1:N)];
m=[(-N+1):(N-1)];
subplot(212)
stem(m,rm);
title('xb自相关');

实验三 IIR数字滤波器的设计
一、 实验目的
1． 掌握双线性变换法及脉冲响应不变法设计IIR数字滤波器的具体设计方法及其原理，熟悉用双线性变换法及脉冲响应不变法设计低通、高通和带通IIR数字滤波器的计算机编程。
2． 观察双线性变换及脉冲响应不变法设计的滤波器的频域特性，了解双线性变换法及脉冲响应不变法设计的特点。
3． 熟悉巴特沃思滤波器、切比雪夫滤波器和椭圆滤波器的频率特性。

二、 实验内容：
（1）

fc=300;
fr=200;
Rb=0.8;
At=20;
T=0.001;
wc=2*1000*tan(2*pi*300/(2*1000));
wt=2*1000*tan(2*pi*200/(2*1000));
[N,wn]=cheb1ord(wc,wt,0.8,20,'s');
[B,A]=cheby1(N,0.8,wn,'high','s');
[num,dem]=bilinear(B,A,1000);
[h,w]=freqz(num,dem);
f=w/pi*500;
plot(f,20*log10(abs(h)));
axis([0,500,-300,10]);
grid;
xlabel('频率 /HZ');
ylabel('幅度 /dB');

（2）
双线性变换法：
Wp=0.2*1000*pi;
Ws=0.3*1000*pi;
Rp=1;
At=25;
Fs=1000;
Ts=1/Fs;
wp=2*tan(Wp/2)/Ts;
ws=2*tan(Ws/2)/Ts;
[N,wn]=buttord(Wp,Ws,Rp,At,'s');
fprintf('滤波器阶数N=%.0f\n',N);
[b,a]=butter(N,wn,'s');
[bz,az]=bilinear(b,a,Fs);
disp('分子系数ba');
fprintf('%.4e',bz);
fprintf('\n');
disp('分母系数aa');
fprintf('%.4e',bz);
fprintf('\n');
w=linspace(0,pi,1024);
h=freqz(bz,az,w);
plot(w/pi,20*log10(abs(h)));
axis([0 1 -350 50]);
grid;
xlabel('频率 ');
ylabel('幅度(dB)');

脉冲响应不变法：
wp=0.2*1000*pi;
ws=0.3*1000*pi;
Rp=1;
At=25;
Fs=1000;
Ts=1/Fs;
Wp=wp*Fs;
Ws=ws*Fs;
[N,Wn]=buttord(Wp,Ws,Rp,At,'s');
[b,a]=butter(N,Wn,'s');
[bz,az]=impinvar(b,a);
w=linspace(0,pi,1024);
h=freqz(bz,az,w);
plot(w/pi,20*log10(abs(h)));
grid;
xlabel('频率');
ylabel('幅度(dB)');

（3）
Wp=1.2*1000*pi;
Ws=2*1000*pi;
Rp=0.5;
At=40;
Fs=8000;
Ts=1/Fs;
wp=2*tan(Wp/2)/Ts;
ws=2*tan(Ws/2)/Ts;
[N,wn]=buttord(Wp,Ws,Rp,At,'s');
fprintf('滤波器阶数N=%.0f\n',N);
[b,a]=butter(N,wn,'s');
[bz,az]=bilinear(b,a,Fs);
w=linspace(0,pi,1024);
h=freqz(bz,az,w);
plot(w/pi,20*log10(abs(h)));
axis([0 1 -350 50]);
grid;
xlabel('频率/kHZ ');
ylabel('幅度(dB)');

巴特沃思滤波器阶数N=12
>>

clear;
Wp=1.2*1000*pi;
Ws=2*1000*pi;
Rp=0.5;
At=40;
Fs=8000;
Ts=1/Fs;
[N,Wn]=cheb1ord(Wp,Ws,Rp,At,'s');
fprintf('滤波器阶数N=%.0f\n',N);
[z0,p0,k0]=cheb1ap(N,Rp);
b0=k0*real(poly(z0));
a0=real(poly(p0));
[H,W]=freqs(b0,a0);
plot(W*Wn,20*log10(abs(H)/max(abs(H)))),grid ;
xlabel('频率/HZ');
ylabel('幅度(dB)');

切比雪夫滤波器阶数N=6
>>
clear;
Wp=1.2*1000*pi;
Ws=2*1000*pi;
Rp=0.5;
At=40;
Fs=8000;
Ts=1/Fs;
[N,Wn]=ellipord(Wp,Ws,Rp,At,'s');
fprintf('滤波器阶数N=%.0f\n',N);
[b,a]=ellip(N,Rp,At,Wn,'s');

omega=[0:200:18000*pi];
h=freqs(b,a,omega);
gain=20*log10(abs(h))
plot(omega/(2*pi),gain);
axis([0 9000 -100 10]);
grid;
xlabel('频率/HZ ');
ylabel('幅度(dB)');

椭圆滤波器阶数N=4
（4）

w1=2*30000*tan(2*pi*2000/(2*30000));
w2=2*30000*tan(2*pi*3000/(2*30000));
wr=2*30000*tan(2*pi*6000/(2*30000));
Rp=3;
At=5;
Fs=30000;
[N,wn]=buttord([w1 w2],[1 wr],3,5,'s');

[B,A]=butter(N,wn,'s');
[num1,den1]=impinvar(B,A,Fs);
[h1,w]=freqz(num1,den1);

[B,A]=butter(N,wn,'s');
[num2,den2]=bilinear(B,A,Fs);
[h2,w]=freqz(