信号分析与处理实验报告2(完成)
信号分析与处理实验报告
实验一1、基本信号的表示及可视化(1)单位冲激信号 (t)程序：t=-1:0.001:1; %定义时间向量for i=1:3; %采用循环语句观察i取不同值时的图形dt=1/(i^4);X=(1/dt)*((t>=(-1/2*dt))-(t>=(1/2*dt)));%计算函数值subplot(1,3,i);%将图像分成三部分同时观察stairs(t,X);title('单位冲激信号δ(t)');end（2）单位阶跃信号程序：t=-0.5:0.001:1;%定义时间变量，间隔为0.001S=stepfun(t,0);%定义单位阶跃信号S1=stepfun(t,0.5);%定义单位阶跃延迟信号figure(1);plot(t,S);axis([-0.5 1 -0.2 1.2]);title('单位阶跃信号')%画出图形figure(2);plot(t,S1);axis([-0.5 1 -0.2 1.2]);title('单位阶跃延迟信号')（3）抽样信号f=sin(t)./t;程序：t=-10:0.6:10; %向量t时间范围t=t1:p:t2,p为时间间隔f=sin(t)./t;plot(t,f,'o'); %显示该信号的时域波形title('f(t)=Sa(t)(时间间隔为0.6s)'); %标题xlabel('t') %横坐标标题axis([-10,10,-0.4,1.1]) %横坐标和纵坐标范围（4）单位样值序列和单位阶跃序列A．单位序列δ(k)B．单位阶跃序列ε(k)程序：n1=-10;n2=10;%输入序列的起始点n=n1:n2;k=length(n);x1=zeros(1,k);x1(1,-n1+1)=1;%产生单位样值序列subplot(2,1,1);%绘图stem(n,x1,'filled');x2=ones(1,k);x2(1,1:-n1)=0;subplot(2,1,2);stem(n,x2,'filled');2、信号的频域分析已知周期方波信号0||2()0||22E t f t T t ττ⎧<⎪⎪=⎨⎪<<⎪⎩，当02T τ=， 04T τ=，08T τ=时，画出其幅度谱和相位谱，观察不同周期下，()f t 的频谱图有何区别。
数字信号处理 实验报告 实验二 应用快速傅立叶变换对信号进行频谱分析
数字信号处理实验报告实验二应用快速傅立叶变换对信号进行频谱分析2011年12月7日一、实验目的1、通过本实验，进一步加深对DFT 算法原理合基本性质的理解，熟悉FFT 算法 原理和FFT 子程序的应用。
2、掌握应用FFT 对信号进行频谱分析的方法。
3、通过本实验进一步掌握频域采样定理。
4、了解应用FFT 进行信号频谱分析过程中可能出现的问题，以便在实际中正确应用FFT 。
二、实验原理与方法1、一个连续时间信号)(t x a 的频谱可以用它的傅立叶变换表示()()j t a a X j x t e dt +∞-Ω-∞Ω=⎰2、对信号进行理想采样，得到采样序列()()a x n x nT =3、以T 为采样周期，对)(n x 进行Z 变换()()n X z x n z +∞--∞=∑4、当ωj ez =时，得到序列傅立叶变换SFT()()j j n X e x n e ωω+∞--∞=∑5、ω为数字角频率sT F ωΩ=Ω=6、已经知道：12()[()]j a m X e X j T T Tωωπ+∞-∞=-∑ ( 2-6 )7、序列的频谱是原模拟信号的周期延拓，即可以通过分析序列的频谱，得到相应连续信号的频谱。
（信号为有限带宽，采样满足Nyquist 定理）8、无线长序列可以用有限长序列来逼近，对于有限长序列可以使用离散傅立叶变换（DFT ）。
可以很好的反映序列的频域特性，且易于快速算法在计算机上实现。
当序列()x n 的长度为N 时，它的离散傅里叶变换为：1()[()]()N knN n X k DFT x n x n W-===∑ 其中2jNN W eπ-=，它的反变换定义为：101()[()]()N knN k x n IDFT X k X k W N --===∑比较Z 变换式 ( 2-3 ) 和DFT 式 ( 2-7 )，令kN z W -=则1()()[()]|kNN nkN N Z W X z x n W DFT x n ---====∑ 因此有()()|kNz W X k X z -==k N W -是Z 平面单位圆上幅角为2kNπω=的点，也即是将单位圆N 等分后的第k 点。
信号分析与处理实验报告
时域采样定理实验日志1实验题目：时域采样定理实验目的：1. 学习掌握matlab 的编程知识及其matalab 在数字信号处理方面常用的12个函数2. 熟悉理想采样的性质，了解信号采样前后的频谱变化，加深对采样定理的理解。
实验要求：1;能对时域上的信号转化为频域上的信号。
2:学会用两种不同的方法保存matlab程序。
3：对给定的模拟信号X a(t) = Ae sin(Ω0t )U (t ) 进行采样！（fm=500），然后把fm改为1500.看看图像的变化实验主要步骤：1：打开matlab编程软件2：输入代码>>n = 0:50-1;>>fs =1000;>>string ='1000';>>Xa=444.128*ex p((-222.144)*n/fs).*sin(222.144*n/fs);>>DF T(Xa,50,str i ng);实验结果：思考题：1，观察实验内容1 中，在分别采用500Hz，1000Hz，1500Hz采样后，对所得的到的信号Xa(n) 绘制的3 个幅频特性曲线有何不同，并分析为什么？结合时域采样定理的内容对图形进行解释；心得体会：时域离散系统及其响应实验日志21:实验题目时域离散系统及其响应2实验目的：1. 继续熟悉掌握m atlab 的使用和编程。
2. 熟悉掌握时域离散系统的时域特性。
并验证时域卷积定理3：实验要求1：能熟练的运用DFT函数，以及DFT函数调用的返回值，并且能运用conv卷积函数2:对于不同的定义域之间的函数卷积，取相同的定义域部分4：实验步骤编写DFT代码编写conv卷积代码：5：实验结果系统h1 (n) = δ(n) + 2.5δ(n −1) + 2.5δ(n −2) + δ(n −3)的图像及其作DFT变换的图像如下输入为x1 (n) = δ(n)的图像及其DFT变换的图像是关于h1和x1作卷积后的变换，以及DFT变换后的图像输入信号为x2 (n) = R10 (n)上述信号作卷积及其DFT变换后的图像输入信号为x3 (t ) = R5 (n)图像及其DFT变换6：实验思考1．比较y1 (n) 和h1 (n) 的时域和频域特性，注意它们之间有无差别，用所学理论解释所得结果。
统计信号分析与处理实验报告
实验2 随机过程的计算机模拟一、实验目的1、给定功率谱（相关函数）和概率分布，通过计算机模拟分析产生相应的随机过程；2、通过该随即过程的实际功率谱（相关函数）和概率分布验证该实验的有效性；3、学会运用Matlab 函数对随机过程进行模拟。
二、实验原理1、标准正态分布随机序列的产生方法：利用随机变量函数变换的方法。
设r1，r2为两个相互独立的（0，1）均匀分布的随机数，如果要产生服从均值为m,方差为正态分布的随机数x，则可以按如下变换关系产生：2、正态随机矢量的模拟：设有一 N 维正态随机矢量，其概率密度为为协方差矩阵，且是正定的。
3、具有有理谱的正态随机序列的模拟根据随机过程通过线性系统的理论，白噪声通过线性系统后，输出是正态的，且输出功率谱只与系统的传递函数有关。
利用这一性质，我们可以产生正态随机过程。
如上图所示，输入W(n)为白噪声，假定功率谱密度为G (z) = 1 W ，通过离散线性系统后，输出X (n)是正态随机序列，由于要求模拟的随机序列具有有理谱，则G (z) X 可表示为：其中，G (z) X+ 表示有理谱部分，即所有的零极点在单位圆之内，G (z) X? 表示非有理谱部分，即所有零极点在单位圆之外。
4、满足一定相关函数的平稳正态随机过程的模拟，当已知平稳随机过程的相关函数而要确定该随机过程的模拟算法。
很显然，只要我们设计一个合适的滤波器，使得该白噪声通过滤波器后，输出的功率谱满足上述相关函数的傅里叶变换，就可以模拟得到该随机过程。
三、实验内容1、产生两组相互独立的（0，1）均匀分布的随机数（随机数个数：500）程序及图形如下:clear;x=randn(1,500);y= randn(1,500);subplot(2,1,1);plot(x);title('第二组');subplot(2,1,2);plot(y);title('第一组')2、按照实验原理中的方法产生一组均值为1，方差为1 的正态分布的随机序列（序列长度：500）程序及图形如下：clear;y=1+sqrt(1)*randn(1,500);plot(y);title(‘正态分布，均值方差都为1’)3、画出功率谱密度为G(w)=1/(1.25+cosw) 的功率谱图（一个周期内），采用均匀采样方法，采样点数为500程序及图形如下：clear;w=rand(1,500);M=1.25+cos(w);N=1;G=N./M;plot(G);title('均匀采样功率频谱');5、模拟产生一个功率谱为G(w)=1/(1.25+cosw) 的正态随机序列程序及图形如下:clear;w=randn(1,500);M=1.25+cos(w);N=1;G=N./M;plot(G);title('均匀采样功率频谱');4、实验中所遇到问题及解决方法问题1、对Matlab软件很生疏、编程也不熟悉。
信号分析与处理实验报告
华北电力大学实验报告||实验名称FFT的软件实现实验（Matlab）IIR数字滤波器的设计课程名称信号分析与处理||专业班级：电气化1308 学生姓名：袁拉麻加学号： 2 成绩：指导教师：杨光实验日期： 2015-12-17快速傅里叶变换实验一、实验目的及要求通过编写程序，深入理解快速傅里叶变换算法（FFT）的含义，完成FFT和IFFT算法的软件实现。
二、实验内容利用时间抽取算法，编写基2点的快速傅立叶变换（FFT）程序；并在FFT程序基础上编写快速傅里叶反变换（IFFT）的程序。
三：实验要求1、FFT和IFFT子程序相对独立、具有一般性，并加详细注释；2、验证例6-4，并能得到正确结果。
3、理解应用离散傅里叶变换（DFT）分析连续时间信号频谱的数学物理基础。
四、实验原理：a.算法原理1、程序输入序列的元素数目必须为2的整数次幂，即N=2M,整个运算需要M 级蝶形运算；2、输入序列应该按二进制的码位倒置排列，输出序列按自然序列排列；3、每个蝶形运算的输出数据军官占用其他输入数据的存储单元，实现“即位运算”；4、每一级包括N/2个基本蝶形运算，共有M*N/2个基本蝶形运算；5、第L级中有N/2L个群，群与群的间隔为2L。
6、处于同一级的各个群的系数W分布相同，第L级的群中有2L-1个系数；7、处于第L级的群的系数是(p=1,2,3，…….,2L-1)而对于第L级的蝶形运算，两个输入数据的间隔为2L-1。
b.码位倒置程序流程图开始检测A序列长度nk=0j=1x1(j)=bitget(k,j);j=j+1Yj<m?Nx1=num2str(x1);y(k+1)=bin2dec(x1);clear x1k=k+1c.蝶形运算程序流程图五、程序代码与实验结果a.FFT程序：%%clear all;close all;clc;%输入数据%A=input('输入x(n)序列','s');A=str2num(A);% A=[1,2,-1,4]; %测试数据%%%%校验序列,%n=length(A);m=log2(n);if (fix(m)~=m)disp('输入序列长度错误，请重新输入!');A=input('输入x(n)序列','s');A=str2num(A);elsedisp('输入正确，请运行下一步')end%%%码位倒置%for k=0:n-1for j=1:m %取M位的二进制数%x1(j)=bitget(k,j); %倒取出二进制数%endx1=num2str(x1); %将数字序列转化为字符串%y(k+1)=bin2dec(x1); %二进制序列转化为十进制数%clear x1endfor k=1:nB(k)=A(y(k)+1); %时间抽取序列%endclear A%%%计算%for L=1:m %分解为M级进行运算%LE=2^L; %第L级群间隔为2^L%LE1=2^(L-1); %第L级中共有2^(L-1)个Wn乘数，进行运算蝶运算的两数序号相隔LE1%W=1;W1=exp(-1i*pi/LE1);for R=1:LE1 %针对第R个Wn系数进行一轮蝶运算，共进行LE1次%for P=R:LE:n %每个蝶的大小为LE% Q=P+LE1;T=B(Q)*W;B(Q)=B(P)-T;B(P)=B(P)+T;endW=W*W1;endendB %输出X(k)%%%验证结果：例6-4b.IFFT程序：%%clear all;close all;clc;%输入数据%A=input('输入X(k)序列','s');A=str2num(A);% A=[6,2+2i,-6,2-2i]; %测试数据%%%%校验序列,%n=length(A);m=log2(n);if (fix(m)~=m)disp('输入序列长度错误，请重新输入!');A=input('输入x(n)序列','s');A=str2num(A);elsedisp('输入正确，请运行下一步')end%%%码位倒置%for k=0:n-1for j=1:m %取M位的二进制数%x1(j)=bitget(k,j); %倒取出二进制数%endx1=num2str(x1); %将数字序列转化为字符串%y(k+1)=bin2dec(x1); %二进制序列转化为十进制数%clear x1endfor k=1:nB(k)=A(y(k)+1); %时间抽取序列%endclear A%%%计算%for L=1:m %分解为M级进行运算%LE=2^L; %第L级群间隔为2^L%LE1=2^(L-1); %第L级中共有2^(L-1)个Wn乘数，进行运算蝶运算的两数序号相隔LE1%W=1;W1=exp(-1i*pi/LE1);for R=1:LE1 %针对第R个Wn系数进行一轮蝶运算，共进行LE1次%for P=R:LE:n %每个蝶的大小为LE%Q=P+LE1;T=B(Q)*W;B(Q)=B(P)-T;B(P)=B(P)+T;endW=W*W1;endendB=conj(B); %取共轭%B=B/n %输出x(n)%验证结果：六、实验心得与结论本次实验借助于Matlab软件，我避开了用C平台进行复杂的复数运算，在一定程度上简化了程序，并添加了简单的检错代码，码位倒置我通过查阅资料，使用了一些函数，涉及到十-二进制转换，数字-文本转换，二-文本转换，相对较复杂，蝶运算我参考了书上了流程图，做些许改动就能直接实现。
统计信号分析与处理实验报告
一、实验目的和要求了解随机过程特征估计的基本概念和方法，学会运用Matlab 函数对随机过程进行特征估计,并且通过实验了解不同估计方法所估计出来结果之间的差异。
二、实验环境、内容和方法内容：实验原理设随机序列X(n)、Y(n)为各态历经过程，样本分别为x(n)、y(n)（n=0,1,....N-1）。
1、均值的估计2、方差的估计方差估计有两种情况，如果均值m x 已知，则如果均值未知，那么3、相关函数估计4、功率谱估计功率谱的估计有几种方法，（1）自相关法先求相关函数的估计然后对估计的相关函数做傅立叶变换，（2）周期图法先对序列x(n)做傅立叶变换，则功率谱估计为周期图法是一种非参数谱估计方法，另外还有一种修正的周期图方法，也叫Welch 法，MATLAB 有周期图和Welch 法的谱估计函数。
（3）现代谱估计技术现代谱估计主要有参数谱估计和子空间谱估计。
参数谱估计法是假定待估计功率谱的信号是白噪声驱动线性系统的输出，常用的方法有基于最大墒估计的伯格算法和Yuler-Walk自回归(AR)方法，这些方法是估计线性系统的参数，通常会得到比经典谱估计方法更好的估计。
子空间法也称为高分辨率谱估计或超分辨率谱估计，常用的方法有MUSIC 法和特征矢量法，这些方法是根据相关矩阵的特征分析或特征分解得到对信号的频率分量的估计，特别适合于线谱（即正弦信号）的估计，是低信噪比环境下检测正弦信号的有效方法。
三、实验过程描述第一题1）、产生一组均值为1，方差为4 的正态分布的随机序列（1000 个样本），估计该序列的均值与方差。
2）画出图像的代码如下：n=1:1:1000;w(n)=1+randn(1,1000)*sqrt(4);plot(n,w(n));mean(w(n))ans =1.0711var(w(n))ans =4.05493）运行以上代码，得到如下图:第二题1）按如下模型产生一组随机序列：x(n)=0.8x(n-1)+w(n)，其中w(n)为均值为1，方差为4 的正态分布白噪声序列。
信号分析与处理实验报告
实验一图像信号频谱分析及滤波一：实验原理FFT不是一种新的变化，而是DFT的快速算法。
快速傅里叶变换能减少运算量的根本原因在于它不断地把长序列的离散傅里叶变换变为短序列的离散傅里叶变换，在利用的对称性和周期性使DFT运算中的有些项加以合并，达到减少运算工作量的效果。
为了消除或减弱噪声，提取有用信号，必须进行滤波，能实现滤波功能的系统成为滤波器。
按信号可分为模拟滤波器和数字滤波器两大类。
数字滤波器的关键是如何根据给定的技术指标来得到可以实现的系统函数。
从模拟到数字的转换方法很多，常用的有双线性变换法和冲击响应不变法，本实验主要采用双线性变换法。
双线性变换法是一种由s平面到z平面的映射过程，其变换式定义为：数字域频率与模拟频率之间的关系是非线性关系。
双线性变换的频率标度的非线性失真是可以通过预畸变的方法去补偿的。
变换公式有Ωp=2/T*tan(wp/2)Ωs=2/T*tan(ws/2)二：实验内容1．图像信号的采集和显示选择一副不同彩色图片，利用Windows下的画图工具，设置成200*200像素格式。
然后在Matlab软件平台下，利用相关函数读取数据和显示图像。
要求显示出原始灰度图像、加入噪声信号后的灰度图像、滤波后的灰度图像。
2．图像信号的频谱分析要求分析和画出原始灰度图像、加入噪声信号后灰度图像、滤波后灰度图像信号的频谱特性。
3.数字滤波器设计给出数字低通滤波器性能指标:通带截止频率fp＝10000 Hz，阻带截止频率fs＝15000 Hz，阻带最小衰减Rs＝50 dB，通带最大衰减Rp＝3 dB，采样频率40000Hz。
三：实验程序clear allx=imread('D:\lan.jpg');%原始彩色图像的数据读取x1=rgb2gray(x);%彩色图像值转化为灰度图像值[M,N]=size(x1);%数据x1的长度,用来求矩阵的大小x2=im2double(x1);%unit8转化为double型x3=numel(x2);%计算x2长度figure(1);subplot(1,3,1);imshow(x2);title('原始灰度图')z1=reshape(x2,1,x3);%将二维数据转化成一维数据g=fft(z1);%对图像进行二维傅里叶变换mag=fftshift(abs(g));%fftshift是针对频域的，将FFT的DC分量移到频谱中心K=40000;Fs=40000;dt=1/Fs;n=0:K-1;f1=18000;z=0.1*sin(2*pi*f1*n*dt);x4=z1+z;%加入正弦噪声f=n*Fs/K;y=fft(x4,K);z2=reshape(x4,M,N);%将一维图转换为二维图subplot(1,3,2);imshow(z2);title('加入噪声后')g1=fft(x4);mag1=fftshift(abs(g1));%设计滤波器ws=0.75*pi;wp=0.5*pi;fs=10000;wp1=2*fs*tan(wp/2);ws1=2*fs*tan(ws/2);rs=50;rp=3;% [n,wn]=buttord(wp/pi,ws/pi,rp,rs);% [bz,az]=butter(n,wn);[n,wn]=buttord(wp1,ws1,rp,rs,'s');[z,p,k]=buttap(n);[b,a]=zp2tf(z,p,k);[B,A]=lp2lp(b,a,wn);[bz,az]=bilinear(B,A,fs);[h,w]=freqz(bz,az,128,fs);L=numel(z2);z3=reshape(z2,1,L);x6=filter(bz,az,double(z3));x7=reshape(x6,M,N);subplot(1,3,3);imshow(x7);g2=fft(x6);mag2=fftshift(abs(g2));title('滤波后')%建立频谱图figure(2);subplot(1,3,1);plot(mag);title('原始Magnitude')subplot(1,3,2);plot(mag1);title('加噪声Magnitude')subplot(1,3,3);plot(mag2);title('滤波后Magnitude')figure(3);subplot(1,2,1)plot(w,abs(h));xlabel('f');ylabel('h');title('滤波器幅谱');subplot(1,2,2);plot(w,angle(h));title('滤波器相谱');四：实验结果与分析图一图二分析：由图二可以知道加入噪声后的幅值谱和原始图的幅值谱明显多了两条幅值线，而这两条幅值线就是我们对原始灰度图加入的正弦噪声，而相应的图一中的加噪声后的图与原始图相比，出现了明显的变化。
合工大自动化信号分析与处理实验报告
合肥工业大学电气与自动化工程学院实验报告自动化专业班级学号姓名日期指导教师黄云志共页第页成绩实验一含噪声语音信号频谱分析及滤波实验报告要求：1、实验报告包括四部分：实验原理、实验内容、实验程序、结果分析；分别占实验报告总成绩的10%，10%，40%，40%；2、实验程序及结果分析如内容雷同，均不给分；一、实验原理含噪声语音信号通过低通滤波器，高频的噪声信号会被过滤掉，得到清晰的无噪声语音信号。
二、实验内容录制一段个人自己的语音信号，并对录制的信号进行采样；画出采样后语音信号的时域波形和频谱图；在语音信号中增加正弦噪声信号(自己设置几个频率的正弦信号)，对加入噪声信号后的语音信号进行频谱分析；给定滤波器的性能指标，采用窗函数法和双线性变换设计数字滤波器，并画出滤波器的频率响应；然后用自己设计的滤波器对采集的信号进行滤波，画出滤波后信号的时域波形和频谱，并对滤波前后的信号进行对比试听，分析信号的变化。
给出数字低通滤波器性能指标:如，通带截止频率fp＝10000Hz，阻带截止频率fs＝12000Hz （可根据自己所加入噪声信号的频率进行阻带截止频率设置），阻带最小衰减Rs＝50dB，通带最大衰减Rp＝3dB（也可自己设置），采样频率根据自己语音信号采样频率设定。
三、实验程序1.原始信号fs=22050;%设置采样频率x1=wavread('C:\Users\m\Documents\FFOutput\无标题.wav');%读取声音sound(x1,fs);%播放原始声音y1=fft(x1,length(x1));%FFTf=fs*(0:511)/1024;figure(1);%建立图形plot(x1);%原始信号时域图形title('原始语音信号时域波形')xlabel('时间');ylabel('幅值');axis([0,5*10^4,-2,2]);figure(2);plot(f,abs(y1(1:512)))%做原始语音信号的FFT频谱图xlabel('频率');ylabel('幅度');title('原始语音信号FFT频谱幅频')2.加噪声clear allfs=22050;%设置采样频率x1=wavread('C:\Users\m\Documents\FFOutput\无标题.wav');%读取声音f=fs*(0:511)/1024;t=0:1/22050:(length(x1)-1)/22050;d=[0.5*cos(2*pi*10000*t)]';%加噪声x2=x1+d;sound(x2,22050);%播放加噪音后的信号y2=fft(x2);figure(3)plot(t,x2)%加噪信号时域图形title('加噪声后的时域信号波形'); xlabel('时间');ylabel('幅值');figure(4)subplot(211);plot(abs(x1));title('原始语音信号频谱');xlabel('频率');ylabel('幅值');subplot(212);plot(abs(y2));title('加噪声后的信号频谱'); xlabel('频率');ylabel('幅值');3．设计IIR滤波器clear allfs=22050;wp=0.4*pi;ws=0.7*pi;rp=3;rs=50;x1=wavread('C:\Users\m\Documents\FFOutput\无标题.wav');%读取声音f=fs*(0:22499)/22500;%F=([1:N]-1)*Fs/N;%换算成实际的频率值t=0:1/22050:(length(x1)-1)/22050;%定义噪声信号d=[0.7*sin(3*pi*10000*t)]';x2=x1+d;%加噪声[n,wn]=buttord(wp/pi,ws/pi,rp,rs);%设计IIR滤波器[bz,az]=butter(n,wn/pi);[h,w]=freqz(bz,az,512,fs);%滤波器幅频响应x3=filter(bz,az,x2);%滤波y0=fft(x2,22500);y1=fft(x3,22500);sound(x3,22500);%播放滤波后的声音figure(1)plot(abs(h));%滤波器幅频响应图title('滤波器幅频响应图')figure(2)subplot(2,1,1)plot(t,x2)%画出滤波前的时域图title('滤波前的时域波形');subplot(2,1,2)plot(t,x3);title('滤波后的时域波形');figure(3)subplot(211)plot(f,abs(y0(1:22500)));title('滤波前频谱')subplot(212)plot(f,abs(y1(1:22500)));%滤波后频谱title('滤波后频谱')4.设计FIR滤波器clear allfs=22050;%设置采样频率x1=wavread('C:\Users\m\Documents\FFOutput\无标题.wav');%读取声音f=fs*(0:22499)/22500;%F=([1:N]-1)*Fs/N;%换算成实际的频率值t=0:1/44100:(length(x1)-1)/44100;%定义噪声信号d=[0.5*sin(2*pi*80000*t)]';x2=x1+d;%加噪声N=512;b1=fir1(N,0.3,hanning(N+1));%设计FIR滤波器[h,w]=freqz(b1,1,512,fs);%滤波器幅频响应x3=filter(b1,1,x2);%滤波y0=fft(x2,22500);y1=fft(x3,22500);sound(x3,22500);%播放滤波后的声音figure(1)plot(w,abs(h));%滤波器幅频响应图title('FIR滤波器幅频响应图')figure(2)subplot(211)plot(f,abs(y0(1:22500)));title('滤波前频谱')subplot(212)plot(f,abs(y1(1:22500)));%滤波后频谱title('滤波后频谱')。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
合肥工业大学电气与自动化工程学院 实验报告专业 班级 学号 姓名 日期 指导教师 共 页 第 页实验二 数字滤波器设计实验报告要求：1、实验内容为实验指导书实验三第2题、实验四第1题；2、实验报告包括四部分：实验原理、实验内容、实验程序、结果分析；分别占实验报告总成绩的20%，10%，30%，40%；3、实验程序及结果分析如有内容雷同，均不给分；4、实验程序、结果分析内容可以直接打印。
【 项目一 双线性变换法设计巴特沃斯滤波器 】实验内容：给定待设计的数字高通和带通滤波器的技术指标如下：(1) HP ：Hz f p 400=，Hz f s 300=，Hz F s 1000=，dB p 3=α，dB s 35=α。
(2) BP ：Hz f sl 200=，Hz f 3001=，Hz f 4002=，Hz f sh 500=，Hz F s 2000=，dB p 3=α，dB s 40=α。
试用双线性变换分别设计满足上述要求的巴特沃斯滤波器，给出其系统函数、对数幅频及相频曲线。
实验原理：为了克服脉冲响应不变法中产生频率响应的混叠失真，采用非线性频率压缩方法，将整个频率轴上的频率范围压缩到-π/T ～π/T 之间，再用z =e sT转换到Z 平面上。
也就是说，第一步先将整个S 平面压缩映射到S 1平面的-π/T ～π/T 一条横带里；第二步再通过标准变换关系z =e s 1T将此横带变换到整个Z 平面上去。
这样就使S 平面与Z 平面建立了一一对应的单值关系，消除了多值变换性，也就消除了频谱混叠现象。
由下图看出，在零频率附近，模拟角频率Ω与数字频率ω之间的变换关系接近于线性关系；但当Ω进一步增加时，ω增长得越来越慢，最后当Ω→∞时，ω终止在折叠频率ω=π处，因而双线性变换就不会出现由于高频部分超过折叠频率而混淆到低频部分去的现象，从而消除了频率混叠现象。
将S 平面j Ω轴压缩变换到s1平面j Ω轴上的-π/T 到π/T 一段，可以采用以下变换关系：⎪⎭⎫⎝⎛Ω=Ω2tan 1T 这样±∞=Ω变换到Tπ±=Ω1，0=Ω变换到01=Ω，可将上式写成22221111T jT j T jTjeee e j Ω-ΩΩ-Ω+-=Ω令s j =Ω，11s j =Ω，解析延拓到整个s 平面和s1平面，可得22221111T s T s T s T s eee e s --+-=再将1s 平面通过以下变换关系映射到z 平面，即T s e z 1=从而得到s 平面和z 平面的单值映射关系为1111--+-=z z s ， ss z -+=11一般来说，为了使模拟滤波器的某一频率与数字滤波器的任一频率有对应关系，可引入待定常数c ，⎪⎭⎫ ⎝⎛Ω=Ω2tan 1Tc ()Ts Ts e e c s 1111--+-=将Ts ez 1=代入到上式，可得1111--+-=z zcs ， sc s c z -+=在MATLAB 中，双线性Z 变换可以通过bilinear 函数实现，其调用格式为：[Bz ，Az]＝bilinear(B ，A ，Fs)；其中B ，A 为模拟滤波器传递函数G (s )的分子分母多项式的系数向量，而Bz ，Az 为数字滤波器的传递函数H (z )的分子分母多项式的系数向量。
本次滤波器设计的主要方法是先设计低通模拟滤波器，然后转换为高通、带通数字滤波器。
在设计的全过程的各个步骤，matlab 都提供相应的工具箱函数，使得IIR 数字滤波器设计变得简单。
总的来说，我的设计思路主要有以下两种：方案一：从归一化模拟低通原型出发，先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S 域变换到Z 域，而得到所需类型的数字滤波器。
图1 先频率变换再离散方案二：先进行双线性变换，将模拟低通原型滤波器变换成数字低通滤波器；然后在Z 域内经数字频率变换为所需类型的数字滤波器。
图2 先离散再频率变换以上两种方案都可以，我最后选择了第一种方案进行设计，即先在模拟域内经频率变换成为所需类型的模拟滤波器；然后进行双线性变换，由S 域变换到Z 域，而得到所需类型的数字滤波器。
因为高通，带通滤波器的设计方法是先将要设计的滤波器的技术特性指标通过频率转换成模拟低通滤波器的技术指标，再根据这些性能指标设计出低通滤波器的传递函数，为了设计程序的简洁，故选择了方案一。
下面给出matlab 程序：实验程序： 【高通滤波器】fp=400;fs=300;Rp=3;Rs=35; %衰减参数FS=1000;T=1/(FS);wp=2*pi*fp;ws=2*pi*fs; %把数字域特征换成模拟域特征 Wp=wp/(FS);Ws=ws/(FS); %归一化数字频率 wp2=2*tan(Wp/2)/T;ws2=2*tan(Ws/2)/T; %频率预畸变 [N,Wn]=buttord(wp2,ws2,Rp,Rs,'s');[z,p,k]=buttap(N); %Buttord 低通原型[Bap,Aap]=zp2tf(z,p,k); %零极点转换为传递函数 figure(1)freqs(Bap,Aap); %模拟低通频率响应 title('模拟滤波器（低通原型）频率响应')[Bbs,Abs]=lp2hp(Bap,Aap,Wn); %模拟低通变高通 [numz denz]=butter(N,Wp/pi,'high');Gz_HP=tf(numz,denz,T) %求传递函数 figure(2)freqs(Bbs,Abs);title('模拟滤波器（高通）频率响应')[Bbz,Abz]=bilinear(Bbs,Abs,FS); %双线性变换 figure(3)freqz(Bbz,Abz,512,FS); title('数字滤波器频率响应') [hw,w]=freqz(Bbz,Abz,512); figure(4)plot(w/pi,20*log(abs(hw))); gridaxis([0,1,-740,50])title('详细显示幅频响应')xlabel('w/pi');ylabel('幅度（dB ）');【带通滤波器】fs1=200;fs2=500;fp1=300;fp2=400Rp=3;Rs=40; %衰减参数 FS=2000;T=1/(FS);ws1=2*pi*fs1;ws2=2*pi*fs2;wp1=2*pi*fp1;wp2=2*pi*fp2; %把数字域特征换成模拟域特征Ws1=ws1/(FS);Ws2=ws2/(FS);Wp1=wp1/(FS);Wp2=wp2/(FS); %归一化数字频率WS1=2*tan(Ws1/2)/T;WS2=2*tan(Ws2/2)/T;WP1=2*tan(Wp1/2)/T;WP2=2*tan(Wp2/2)/T; %频率预畸变WP=WP2-WP1;WS=sqrt(WS1*WS2);[N,Wn]=buttord(WP,WS,Rp,Rs,'s');[z,p,k]=buttap(N); %Buttord低通原型[Bap,Aap]=zp2tf(z,p,k); %零极点转换为传递函数figure(1)freqs(Bap,Aap); %模拟低通频率响应title('模拟滤波器（低通原型）频响')[Bbs,Abs]=lp2bp(Bap,Aap,WS,WP); %模拟低通变带通[numz denz]=butter(N,[Wp1/pi Wp2/pi]);Gz_BP=tf(numz,denz,T) %求传递函数figure(2)freqs(Bbs,Abs);title('模拟滤波器（带通）频响')[Bbz,Abz]=bilinear(Bbs,Abs,FS); %双线性变换figure(3)freqz(Bbz,Abz,512,FS);title('数字滤波器频响')[hw,w]=freqz(Bbz,Abz,512);figure(4)plot(w/pi,20*log(abs(hw)));gridaxis([0,1,-740,50])title('详细显示幅频响应')xlabel('w/pi');ylabel('幅度（dB）');实验图表：【高通滤波器】【带通滤波器】说明：为了检验所设计的滤波器的性能，我选择了一个信号函数：X=sin(2*pi*f1*t)+cos(2*pi*f2*t)+sin(2*pi*1*t)选择合适的f1和f2,对输入信号进行了滤波，这时候需要借助Matlab中的y=filter(bz,az,x)函数，附加的程序如下：【高通部分】n=0:180;dt=1/FS;t=n*dt; %时间序列f1=100;f2=550; %输入信号频率x=sin(2*pi*f1*t)+cos(2*pi*f2*t)+sin(2*pi*t); %输入信号figure(5);plot(t,x);axis([0,0.05,-2,2]);grid; %绘制输入信号title('输入信号');xlabel('时间t');ylabel('幅度');y=filter(Baz,Aaz,x); %对输入信号进行低通滤波figure(6);plot(t,y);axis([0,0.05,-2,2]);grid; %绘制输出信号title('低通滤波输出信号');xlabel('时间t');ylabel('幅度');y=filter(Bbz,Abz,x); %对输入信号进行高通滤波figure(7);plot(t,y);axis([0,0.05,-2,2]);grid; %绘制输出信号title('高通通滤波输出信号');xlabel('时间t');ylabel('幅度');【带通部分】y=filter(Bbz,Abz,x); %对输入信号进行带通滤波figure(8);plot(t,y);axis([0,0.1,-2,2]);grid; %绘制输出信号title('带通滤波输出信号');xlabel('时间t');ylabel('幅度');下面附上信号X分别通过高通和带通滤波器的仿真图:对应的传递函数如右图：结果分析：由仿真结果可知，设计出的数字高通，带通的幅频响应均符合各自的特性，从模拟滤波器特性曲线可以看出，所设计的高通滤波器在大于400Hz为通带，其衰减均小于3dB；小于300Hz为阻带，其衰减大于35dB。