matlab信号处理相关函数程序

1、序列的傅里叶变换及其逆变换定义： 其幅度特性为，在Matlab中采用abs函数；相位特性为，在Matlab中采用angle函数。为了方便，考虑在两个周期，例如中2M+1个均匀频率点上计算FT，并且观察其周期性和对称性。为此给出function文件如下，求解FT变换：functionX,w=ft(x,n,k)%X：序列x(n)的傅里叶变换%w：X的自变量%x：要进行傅里叶变换的序列x(n)%n：序列x(n)的位置向量%k：求和区间w=(pi/abs(max(k)/2)*k;X=x*(exp(-1i*pi/abs(max(k)/2).(n*k);使用方法如下：n=-5:5;%序列区间x=(-0.9

2、).n;%序列表达式k=-200:200;%求和区间Xw,w=ft(x,n,k);%求傅里叶变换magX=abs(Xw);%求幅度angX=angle(Xw);%求相位realX=real(Xw);imagX=imag(Xw);subplot(2,2,1)plot(w/pi,magX)%绘制幅度曲线grid ontitle(幅度曲线)xlabel(omega/pi)ylabel(幅度)xmin=0;xmax=2;set(gca,xlim,xmin,xmax,ylimmode,auto,zlimmode,auto); %xmin xmax为范围subplot(2,2,2)plot(w/pi,an

3、gX/pi) %绘制相位曲线grid ontitle(相位曲线)xlabel(omega/pi)ylabel(相位)% angle(X)/pixmin=0;xmax=2;set(gca,xlim,xmin,xmax,ylimmode,auto,zlimmode,auto); %xmin xmax为范围subplot(2,2,3)plot(w/pi,realX) %绘制实部曲线grid ontitle(实部曲线)xlabel(omega/pi)ylabel(实部)xmin=0;xmax=2;set(gca,xlim,xmin,xmax,ylimmode,auto,zlimmode,auto);

4、%xmin xmax为范围subplot(2,2,4)plot(w/pi,imagX) %绘制虚部曲线grid ontitle(虚部曲线)xlabel(omega/pi)ylabel(虚部)xmin=0;xmax=2;set(gca,xlim,xmin,xmax,ylimmode,auto,zlimmode,auto); %xmin xmax为范围序列的DFT及IDFT定义： 离散傅里叶变换的的性质：（1）DFT的共轭对称性。若，则：， 。（2）实序列DFT的性质。若为实序列，则其离散傅里叶变换为共轭对称，即。（3）实偶序列DFT的性质。若为实偶序列，则其离散傅里叶变换为实偶对称，即。（4）实

5、奇序列DFT的性质。若为实奇序列，则其离散傅里叶变换为纯虚奇对称，即。离散傅立叶变换函数 function Xk,k=dft(xn,N)n=0:1:N-1;k=0: N-1;WN=exp(-1j*2*pi/N);nk=n*k;WNnk=WN.(nk);Xk=xn*WNnk; %采用矩阵相乘的方法magX=abs(Xk); k=(0:length(magX)-1)*N/length(magX);离散傅立叶反变换函数 function xn=idft(Xk,N)n=0:1:N-1;k=0:1:N-1;WN=exp(-1j*2*pi/N);nk=n*k;WNnk=WN.(-nk);xn=(Xk*WN

6、nk)/N;使用方法如下：1、序列的傅里叶变换及离散傅里叶变换计算N=5;n=0:4;x=ones(1,5); %产生矩形序列k=0:999;w=(pi/500)*k;X=x*(exp(-j*pi/500).(n*k); %计算序列的傅立叶变换Xe=abs(X); subplot(3,2,1);stem(n,x);ylabel(x(n); subplot(3,2,2);plot(w/pi,Xe);ylabel(|X(ejw)|); %画出序列的傅立叶变换X=dft(x,N); %进行5点DFTmagX=abs(X); k=(0:length(magX)-1)*N/length(magX);su

7、bplot(3,2,3);stem(n,x);ylabel(x(n); subplot(3,2,4);stem(k,magX); axis(0,5,0,5);ylabel(|X(k)|); 2、序列产生及其DFTN=20;n=0:N-1;x=2*cos(0.25*pi*n)+cos(0.65*pi*n);subplot(2,1,1);stem(n,x);title(N=20时的信号);Y=dft(x,N);k1=0:N-1;w1=2*pi/N*k1;subplot(2,1,2);stem(w1/pi,abs(Y);title(信号的频谱);FFT算法计算矢量xn的离散傅立叶变换一维快速傅立叶变

8、换函数 fftXk=fft(xn,N) xn为时域序列向量，N是DFT变换区间长度。当N大于xn的长度时，fft函数自动在xn后面补零。函数返回xn的N点DFT变换结果向量Xk。当N小于xn的长度时，fft函数计算xn的前面N个元素构成的N长序列的N点DFT，忽略xn后面的元素。当x为矩阵时，Xk为矩阵xn的每一列的FFT ,fft 函数按类似的方法处理列长度。当xn的长度为2的幂次方时，则fft采用基2的FFT算法，否则采用稍慢的混合基算法。一维快速傅立叶反变换 ifftxn=ifft(Xk,N) Xk的N点IFFTxn=ifft(Xk) 用于计算矢量Xk的IFFT使用方法如下：M=26;N

9、=32;n=0:M;xa=0:M/2;xb=ceil(M/2)-1:-1:0;xn=xa,xb;%产生M长三角波序列x(n)Xk=fft(xn,512);%512点FFTX32k=fft(xn,32);%32点FFTx32n=ifft(X32k);%32点IFFT得到x32(n)X16k=X32k(1:2:N);%隔点抽取X32k得到X16kx16n=ifft(X16k,N/2); %16点IFFT得到x16(n)连续信号的傅里叶变换及其逆变换定义： Fw=fourier(ft,t,w)%求时域函数ft的Fourier变换FwFt=ifourier(Fw,w,t)%求频域函数Fw的Fourie

10、r逆变换ft使用方法如下：syms t wut=heaviside(t);%单位阶跃函数Uw=fourier(ut,t,w);%傅里叶变换SUt=simple(Uw);%化简Inv_ut=ifourier(Uw,w,t);%傅里叶逆变换ezplot(w,Uw)%绘制傅里叶变换曲线title(傅里叶变换曲线)grid on因果序列的Z变换及其逆变换定义： Xz=ztrans(xn,n,z)%求时域序列xn的Z变换Xzxn=iztrans(Xz,z,n)%求频域序列Xz的Z逆变换xn使用方法如下：syms n w T zxn=sin(w*n*T);%定义序列Xz=ztrans(xn,n,z);%Z

11、变换Inv_xn=iztrans(Xz,z,n);%Z逆变换连续信号的单边拉普拉斯变换及其逆变换定义：Fs=laplace(ft,t,s)%求时域函数ft的Laplace变换Fsft=ilaplace(Fs,s,t)%求频域函数Fs的Laplace逆变换ft使用方法如下：syms t ssyms a positiveft=dirac(t);%单位冲激函数Fs=laplace(ft,t,s);%拉普拉斯变换Inv_ft=ifourier(Fs,s,t);%拉普拉斯逆变换线性常系数差分方程的递推求解 求差分方程的零状态响应和全响应函数yn=filter(B,A,xn) 计算系统对输入信号向量xn的

12、零状态响应输出信号向量yn，yn与xn长度相等，其中，B和A是差分方程的系数向量，即 ,其中，如果，则filter用 对系数向量B和A归一化。yn=filter(B,A,xn,xi) 计算系统对输入信号向量xn的全响应输出信号yn。其中，xi是等效初始条件的输入序列，所以xi是由初始条件确定的。xi=filtic(B,A,ys,xs) 其中，ys和xs是初始条件向量：ys=y(-1),y(-2),y(-3),y(-N)，xs=x(-1),x(-2),y(-3),x(-M)。如果xn是因果序列，则xs=0，调用时可缺省xs。使用方法如下：%调用filter解差分方程y(n)-ay(n-1)=x(

13、n)a=0.8;ys=1; %设差分方程系数a=0.8，初始状态：y(-1)=1xn=1,zeros(1,30); %x(n)=单位脉冲序列，长度N=31B=1;A=1,-a; %差分方程系数xi=filtic(B,A,ys); %由初始条件计算等效初始条件的输入序列xiyn=filter(B,A,xn,xi); %调用filter解差分方程，求系统输出函数y(n)n=0:length(yn)-1;stem(n,yn,.)xlabel(n);ylabel(y(n)离散系统函数定义： 系统的因果性：其单位脉冲响应h(n)是因果序列，满足n0时，h(n)=0或 其系统函数H(z)的收敛域包含无穷点

14、系统的稳定性：其单位脉冲响应h(n)绝对可和，即 或 其系统函数H(z)的收敛域包含单位圆系统稳定性判定函数 stab(A)function stab(A)%stab：系统稳定性判定函数，A是H(z)的分母多项式系数向量disp(系统极点为：)P=roots(A)%求H(z)的极点，并显示disp(系统极点模的最大值为：)M=max(abs(P)%求所有极点模的最大值，并显示if MN error(N mustbe = the length of x) %要求移位周期大于信号长度endx=x zeros(1,N-length(x);n=0:1:N-1;n=mod(n-m,N);y=x(n+1)

15、; 时域循环卷积定理：设有限长序列为和，它们的N点DFT分别为和，如果，则其IDFT为两序列的循环卷积。计算两序列的循环卷积函数 function y=circonv(x1,x2,N)%循环卷积的长度N应大于等于max(length(x1),length(x2)if length(x1)N error(N must be = the length of x1)endif length(x2)N error(N must be = the length of x2)endx1=x1 zeros(1,N-length(x1); %将序列补0成为N长序列x2=x2 zeros(1,N-length(

16、x2);m=0:1:N-1;x3=x2(mod(-m,N)+1); %该语句的功能相当于序列翻褶，延拓，取主值序列H=zeros(N,N);for n=1:1:N %得到x3序列的循环移位矩阵 H(n,:)=cirshftt(x3,n-1,N); endy=x1*H %用矩阵相乘的方法得到结果信噪比设纯净信号为，噪声信号为，带噪信号为，则信噪比的定义如式（4-2）所示，单位为dB。 （4-2）function y=snr(x0,xn)%计算信噪比，单位dB%x0是纯净信号,xn是带噪信号p0=sum(abs(x0).2); %sum函数为求和y=10*log10(p0/(sum(abs(xn-

17、x0).2);序列的基本运算%1.单位取样序列 x(n)=delta(n-n0) 要求n1=n0=n2 function x,n=impseq(n0,n1,n2)n=n1:n2; x=(n-n0)=0;%2.单位阶跃序列 x(n)=u(n-n0) 要求n1=n0=0;%3.信号加 y(n)=x1(n)+x2(n)function y,n=sigadd(x1,n1,x2,n2)%find函数：找出非零元素的索引号%x1:第一个序列的值，n1:序列x1的位置向量%x2:第二个序列的值，n2:序列x2的位置向量n=min(min(n1),min(n2):max(max(n1),max(n2);y1=

18、zeros(1,length(n); y2=y1;y1(find(n=min(n1)&(n=min(n2)&(n=min(n1)&(n=min(n2)&(n=max(n2)=1)=x2;y=y1.*y2;%5.移位 y(n)=x(n-n0)function y,n=sigshift(x,m,n0)%m：序列x的位置向量n=m+n0; y=x;%6.翻褶 y(n)=x(-n)function y,n=sigfold(x,n)y=fliplr(x); n=-fliplr(n);特殊的连续函数单位冲激函数 dirac(t)使用方法如下：syms t ndirac(t)delta=sym(kronec

19、kerDelta(n,0);%定义单位冲激单位阶跃函数 heaviside(t)使用方法如下：syms theaviside(t)录音函数function record(n,fs,file) %录音函数%设置n录音时间，采样率fs ,保存文件名filefprintf(Press any key to start %g seconds of recording.n,n);pause;fprintf(Recording.n);y=wavrecord(n*fs,fs,int16); %录音 fprintf(Finished recording.n);wavwrite(y,fs,file); %保存为

20、文件为file fprintf(The recording is saved as );fprintf(file);fprintf(n);读取音频y = wavread(filename) 加载音频文件文件名指定的字符串,返回y的抽样数据。如果不包括一个扩展文件名,wavread默认文件后缀为. wav。y, Fs = wavread(filename) 返回的采样率(Fs)赫兹用于编码中的数据文件。y, Fs, nbits = wavread(filename) 返回每个抽样的比特数(nbits)。写入音频wavwrite(y,filename) 把在变量y中的存储数据写入文件名为filen

21、ame的文件。文件名输入是一个用单引号括起来的字符串。默认数据采样率为8000 Hz和被认为是16位。每一列代表一个单独的数据通道。因此,立体数据应该指定为一个矩阵有两个列。wavwrite(y,Fs,filename) 把在变量y中的存储数据写入文件名为filename的文件。数据的采样率为Fs赫兹和被认为是16位。wavwrite(y,Fs,N,filename) 把在变量y中的存储数据写入文件名为filename的文件。数据采样率为Fs赫兹和N位,其中N是8,16、24、32。播放音频sound(y,Fs) 以采样率Fs播放音频信号y。如果你不指定采样率,默认为8192赫兹。单通道(单声

22、道)音频,y是一个米长、1列向量,m是音频样本的数量。如果你的系统支持立体声播放,y是一个m-by-2矩阵,第一列对应于左声道,和第二列对应于右声道。声音功能假设y包含-1和1之间的浮点数,和剪辑范围之外的值。sound(y,Fs,bits) 指定的位深度(即精度)的样本值。位深度的可能的值取决于可用的音频硬件系统上。大多数平台支持8位或16位的深度。如果你不指定位,声音功能在一个8位深度。wavplay(y,Fs) 播放存储的音频信号向量y。Fs是整数采样率，单位是Hz。Fs的默认值是11025赫兹。wavplay只支持1 或2声道音频(单声道或立体声)信号。在立体声,y必须是一个两列矩阵。

23、使用方法如下：y,fs=wavread(sound.wav);%读取音频文件L=length(y);%采样点数T=L/fs;%采样时间fprintf(按任意键播放原录音及波形n);pause;plot(y); %画图 title(原信号8000Hz波形)wavplay(y,fs);fprintf(录音频率：%g Hzn,fs);fs1=16000;L1=T*fs1;B1=L/L1;y1=0*(1:L1);for i=1:L1 y1(i)=y(ceil(i*B1);%重新采样endwavwrite(y1,fs1,sound1.wav); %保存为文件为sound1.wav y1=wavread(

24、sound1.wav);%读取音频文件fprintf(按任意键播放采样频率 %g Hz录音及波形n,fs1);pause;plot(y1); %画图 title(采样信号16000Hz波形)wavplay(y1,fs1); fs2=8000;L2=T*fs2;B2=L/L2;y2=0*(1:L2);for i=1:L2 y2(i)=y(ceil(i*B2);%重新采样endwavwrite(y2,fs2,sound2.wav); %保存为文件为sound2.wav y2=wavread(sound2.wav);%读取音频文件fprintf(按任意键播放采样频率 %g Hz录音及波形n,fs2);pause;plot(y2); %画图 title(采样信号8000Hz波形)wavplay(y2,fs2);示例1、令，绘制其傅立叶变换。用不同频率对其进行采样，分别画出。Dt=0.00005; %步长为0.00005st=-0.005:Dt:0.005; xa=exp(-1000*abs(t);%取时间从-0.005s到0.005s这段模拟信号Wmax=2*pi*2000; %信号最高频率为2 *2000K=500; %频域正半轴取500个点进行计