实验四信号与系统(数字信号部分2)

实验四信号的频域分析实验目的：1、掌握FFT和IFFT；2、掌握功率密度普；3、掌握数字滤波；一、实验讲义信号的频域分析在实际中有着广泛的应用，是LTI系统分析的基础。对于连续信号的分析，如周期信号的傅里叶级数、非周期信号的傅里叶变换等，都需要利用其相应的公式进行分解或变换，以便获得傅里叶级数和傅里叶变换的表达形式，然后利用绘图工具绘制其频谱。对于连续信号而言，MATLAB无法提供专门的函数进行频谱分析和处理。因此，本节只讨论离散信号和离散系统的频谱分析，并介绍频域分析的应用（以数字滤波为例）。1离散傅里叶变换及其逆变换傅里叶变换可以将时域信号转换为频域信号，以便分析信号的频域特性，其逆变换则把频域信号转换为时域信号。傅里叶变换的原理是把一个时域信号分解成用不同频率的正弦信号（或复信号）线性组合的形式，这样时域信号所包含的频率成分就一目了然。离散傅里叶变换DFT是数字信号分析的主要工具，快速傅里叶变换是对离散傅里叶变换DFT进行快速计算的有效算法。在MATLAB中，利用函数fft和ifft分别计算一维信号的离散傅里叶变换和其逆变换。对于二维信号，使用的离散傅里叶变换和离散傅里叶变换逆变换函数则为fft2和ifft2。函数fft和ifft的调用格式为：y=fft(x,n)；Y=ifft(X,n)说明：（1）x和X分别为待变换的输入向量，x为时域信号，X为频域信号。（2）n表示进行变换的点数，可以缺省，若输入的序列比n短，则fft和ifft用0填充序列，使其长度为n；若输入的序列比n长，则截短输入序列。（3）n缺省时，变换的点数为输入序列的长度。（4）n为2的幂数时，计算速度最快。例1对信号()sin100sin240xttt进行傅里叶变换，然后对变换后的序列进行傅里叶逆变换，并绘制它们的图像。N=512;%N为采样点数T=1;%T为采样时间终点t=linspace(0,T,N);%给出N个采样时间ti(I=1:N)x=sin(2*pi*50*t)+5*cos(2*pi*120*t);%求各采样点样本值xfigure(1),plot(x(1:30));%绘制输入信号的部分波形dt=t(2)-t(1);%采样周期f=1/dt;%采样频率(Hz)X=fft(x);%计算x的快速傅立叶变换XF=X(1:N/2+1);%F(k)=X(k)(k=1:N/2+1)f=f*(0:N/2)/N;%使频率轴f从零开始figure(2),plot(f,abs(F),'-*')%绘制振幅-频率图xlabel('Frequency');ylabel('|F(k)|');z=ifft(X);%傅里叶逆变换Z=real(z);%逆变换之后的信号是复信号，去其实部即可figure(3),plot(Z(1:30));%绘制逆变换之后的信号程序运行结果如图1、图2和图3所示。051015202530-6-4-20246图1输入信号的时域图050100150200250300020040060080010001200Frequency|F(k)|图2输入信号的频谱051015202530-6-4-20246图3逆傅里叶变换后的信号可以看出，1和2是完全一致的，说明傅里叶变换和傅里叶逆变换是严格对应的。从图2可知，频谱在50Hz和120Hz有两个冲激，说明原时域信号中同时含有频率为50Hz和120Hz的两个周期信号。2信号的功率密度谱信号的功率密度谱是信号功率谱密度值在频率范围内的分布，具体表现就是不同频段上波形幅值的变化。功率谱是进行频率分析的基础，在许多工程方面有着重要的意义。在MATLAB中，求取信号的功率密度谱用函数psd来实现，该函数的调用格式为：[Pxx,f]=psd(Xn,nfft,Fs,window,noverlap)说明：（1）Pxx为输入信号Xn的功率谱密度数值序列，f为与Pxx对应的频率序列。（2）Xn为输入的时域信号，nfft为计算FFT的单位宽度，FS为采样频率，window为声明窗函数的类型，noverlap是处理Xn混叠的点数。（3）使用psd函数计算信号的功率密度谱的基本原理是：把输入时域信号Xn分成许多连续的区域，对每个区域加窗（窗函数的类型由window设置），然后做FFT计算，每两个相邻的区域之间有个重叠的区域（大小由noverlap设置）。把每个区域计算所得的序列先做模计算，然后对所有区域的序列进行求和，最后用这个和除以计算过的功率分布范围就获得到了信号的功率密度谱。例2在时域信号()sin(120)2sin(320)xttt掺入随机噪声，并绘制该信号的功率密度谱。t=0:0.001:1;x=sin(2*pi*60*t)+2*sin(2*pi*160*t)+randn(size(t));%信号加噪声nfft=256;%设置nfft值Fs=1000;%设置采样频率window=hanning(nfft);%设置窗函数为汉宁窗，窗宽为256noverlap=128;%混叠宽度为128点[Pxx,f]=psd(x,nfft,Fs,window,noverlap);%计算功率密度谱plot(t(1:200),x(1:200));%绘制输入信号figure,plot(f,Pxx)%绘制功率密度谱程序运行结果如图4所示。050100150200250300350400450500020406080100120140160180图4信号的功率密度谱从图4可以看出，经过psd函数的计算，随机噪声的能量被大大地弱化，而周期信号的能量得到了加强。3信号的互相关功率密度谱信号的相关分为自相关和互相关两种类型，分别说明一个信号与自己或另外一个信号之间在频域上的相似性。信号的互相关功率密度谱（如果两个信号完全相同，则为自相关功率密度谱）在故障诊断和状态预测等方面有广泛的应用，例如检测振动信号中是否有周期信号、检测零部件裂缝的位置及孔洞的大小等。在MATLAB中，求取信号的互相关功率密度谱用函数csd来实现，该函数的调用格式为：[Pxy,f]=csd(x,y,nfft,Fs,window,noverlap)说明：（1）x和y分别为两个时域的输入信号，nfft为计算FFT的单位宽度，Fs为采样频率，window是指窗函数的类型，noverlap是处理Xn混叠的点数。（2）Pxy为互相关功率密度谱的数值序列，f为相应的频率轴。（3）该函数与功率密度谱函数psd的原理基本相似。（4）该函数的主要作用是用来辨认信号中的周期成分。例3求一个含有频率为100Hz正弦波的随机噪声信号的自相关功率密度谱，并绘制其图像。t=0:0.001:1;x=sin(2*pi*100*t)+randn(size(t));%含有频率为100Hz正弦波的随机信号nfft=256;%设置nfft值Fs=1000;%设置采样频率window=hanning(256);%设置窗函数为汉宁窗，窗宽为256noverlap=128;%混叠宽度为128点[Pxx,f]=csd(x,x,nfft,Fs,window,noverlap);%计算自相关功率密度谱subplot(1,2,1),plot(t(1:200),x(1:200));title(‘含有正弦信号的随机噪声’);%绘制输入信号subplot(1,2,2),plot(f,Pxx);title(‘功率密度谱’)%绘制功率密度谱程序运行结果如图5所示。00.050.10.150.2-4-3-2-101234含有正弦信号的随机噪声0200400600051015202530354045功率密度谱图5含有正弦信号的随机噪声及其自相关功率密度谱由图5可知，自相关功率密度值都大于0，事实上自相关功率谱还是偶函数。从时域图像中很难看出此输入信号中含有周期成分，更不要说确定周期成分的周期，但是从自相关谱中很容易就看出这个杂乱无章的信号中含有频率为100Hz的周期成分。例4试绘制两个时域信号之间的互相关功率谱。t=0:0.001:1;x=sin(2*pi*300*t)+2*sin(2*pi*160*t)+randn(size(t));y=sin(2*pi*100*t)+2*sin(2*pi*180*t)+randn(size(t));nfft=256;%设置nfft值Fs=1000;%设置采样频率window=hanning(256);%设置窗函数为汉宁窗，窗宽为256noverlap=128;%混叠宽度为128点[Pxy,f]=csd(x,y,nfft,Fs,window,noverlap);%计算互相关功率密度谱plot(f,Pxy)%绘制功率密度谱图像程序运行结果如图6所示：050100150200250300350400450500-3-2-101234567图6两个信号的互相关功率谱两个信号分别含有频率为160Hz、300Hz和100Hz、180Hz的正弦波，从信号的互相关功率谱可以明显观察到这些频率的周期信号。从图9-20中可知，互相关功率谱不全为正值。4数字滤波滤波器设计的目的是对信号进行数据序列的频率变换，去掉不需要的频率成分，变换成输出序列。一般有低通（滤除高频成分）、高通（滤除低频成分）、带通（滤除通频带两端的频率成分）和带阻（滤除声明的阻频带部分）等几种滤波器。一般来讲，数字滤波器的输出序列y(n)的Z变换Y(z)与滤波器的输入序列x(n)的Z变换X(z)是相关联的，常常表述如下：11(1)(2)...(1)()()()()(1)(2)...(1)mnbbzbmzYzHzXzXzaazanz（1）其中，H(z)是滤波器的传递函数，常量a(i)和b(i)是滤波器的系数，而滤波器的阶次是m和n中的较大值，滤波器的系数分别以两个行向量a和b的形式被存储起来。在MATLAB中，滤波运算用函数filter来实现，该函数的调用格式为：y=filter(b,a,x)其中，a、b是滤波器的系数向量。滤波器的输出序列y和输入序列x的长度（即采样点数）是相等的。在MATLAB中，不同的滤波器类型是通过不同的函数计算来实现的。滤波器有很多类型，其对应的函数也各不相同，下面介绍其中的一类滤波器即Butterworth滤波器。当计算出Butterworth滤波器的系数后，即可用之来实现信号的滤波。例5在正弦信号()sin50xtt中加入随机噪声，采样频率为1000Hz，请滤除其中频率为30Hz以上的噪声。t=0:0.001:0.5;x=sin(50*pi*t)+randn(size(t));%生成输入序列[b,a]=butter(10,30/500);%计算滤波器系数y=filter(b,a,x);%进行数字滤波subplot(1,2,1),plot(t,x);axis([00.5-44]);title(‘含噪信号’);%绘制滤波前信号的图像subplot(1,2,2),plot(t,y);axis([00.5-22]);title(‘滤波结果’);%绘制滤波后信号的图像程序运行结果如图7所示：00.10.20.30.4-4-3-2-101234含噪信号00.10.20.30.4-2-1.5-1-0.500.511.52滤波结果图7输入序列及其滤波结果程序说明：（1）butter(10,30/500)是一个滤波器函数，它的第一个参数表示这个滤波器的阶数是10，此处的阶数为估计值。在不同的情况下应当选取什么样的阶次，MATLAB里面有专门的函数来计算这个数值。butter函数的阶次估计函数为buttord，其具体的用法此处不予介绍。（2）butter(10,30/500)的第二个参数是滤波器的标准化频率。在MATLAB中，包括butter在内的所有滤波器函数的截止频率都采用标准化频率，即以Hz为单位的截止频率除以采样频率的一半，所以标准频率是一个分布在区间[0,1]的比值。（3）从图7可知，信号中的高频成分被滤除了，但是由于其中仍然含有频率为30Hz以下的随机噪声，滤波后的信号不再是标准的正弦信号。