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=

2、(-0.9).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','

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

4、); %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'); %xmin xmax为范围subplot(2,2,4)plot(w/pi,imagX) %绘制虚部曲线grid

5、 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的性质。若为实

6、偶序列，则其离散傅里叶变换为实偶对称，即。（4）实奇序列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);

7、nk=n'*k;WNnk=WN.(-nk);xn=(Xk*WNnk)/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)

8、; %进行5点DFTmagX=abs(X); k=(0:length(magX)'-1)*N/length(magX);subplot(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

9、);k1=0:N-1;w1=2*pi/N*k1;subplot(2,1,2);stem(w1/pi,abs(Y);title('信号的频谱');FFT算法计算矢量xn的离散傅立叶变换一维快速傅立叶变换函数 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的幂次

10、方时，则fft采用基2的FFT算法，否则采用稍慢的混合基算法。一维快速傅立叶反变换 ifftxn=ifft(Xk,N) Xk的N点IFFTxn=ifft(Xk) 用于计算矢量Xk的IFFT使用方法如下：M=26;N=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);

11、%16点IFFT得到x16(n)连续信号的傅里叶变换及其逆变换定义： Fw=fourier(ft,t,w)%求时域函数ft的Fourier变换FwFt=ifourier(Fw,w,t)%求频域函数Fw的Fourier逆变换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=zt

12、rans(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变换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=la

13、place(ft,t,s);%拉普拉斯变换Inv_ft=ifourier(Fs,s,t);%拉普拉斯逆变换线性常系数差分方程的递推求解 求差分方程的零状态响应和全响应函数yn=filter(B,A,xn) 计算系统对输入信号向量xn的零状态响应输出信号向量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是初始条件向

14、量：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(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

15、)n=0:length(yn)-1;stem(n,yn,'.')xlabel('n');ylabel('y(n)')离散系统函数定义： 系统的因果性：其单位脉冲响应h(n)是因果序列，满足n<0时，h(n)=0或 其系统函数H(z)的收敛域包含无穷点系统的稳定性：其单位脉冲响应h(n)绝对可和，即 或 其系统函数H(z)的收敛域包含单位圆系统稳定性判定函数 stab(A)function stab(A)%stab：系统稳定性判定函数，A是H(z)的分母多项式系数向量disp('系统极点为：')P=roots(A)%求H(z)

16、的极点，并显示disp('系统极点模的最大值为：')M=max(abs(P)%求所有极点模的最大值，并显示if M<1 disp('系统稳定')else disp('系统不稳定')end、 极点及频率响应曲线zplane(B,A) 绘制出系统函数H(z)的零极点图。其中B和A为系统函数H(Z)的分子和分母多项式系数向量。假设 则B=B(1),B(2),.,B(M+1), A=A(1),A(2),.,A(N+1)H=freqz(B,A,w) 计算由向量w指定的数字频率点上数字滤波器H(z)的频率响应H()，结果存于H向量中。H,w=freqz

17、(B,A,M) 计算出M个点上的频率响应，存于H向量中，M个频率响应存于w向量中，频率范围自动均匀设在0,上。H,w=freqz(B,A,M,whole) 频率范围自动均匀设在0,2上。freqz(B,A) 自动选取512个频率点计算，频率范围自动均匀设在0,上。不带输出向量的freqz函数将自动绘出幅频响应和相频响应曲线，频率和相位是线性坐标，幅频响应为对数坐标。使用方法如下：B=1,0,0,0,0,0,0,0,-1;A=1;%设置系统函数系数向量B和Asubplot(2,2,1)zplane(B,A);%绘制零极点图title('零极点图')H,w=freqz(B,A);%

18、计算频率响应subplot(2,2,2)plot(w/pi,abs(H)/max(abs(H)%绘制幅频响应曲线xlabel('omega/pi')ylabel('幅度')title('幅频特性曲线')axis(0,1,0,2.5)grid onsubplot(2,2,4)plot(w/pi,angle(H)%绘制相频响应曲线xlabel('omega/pi')ylabel('相位')title('相频特性曲线')grid on实序列分解函数 % 实序列可分解为共扼对称分量 % 和共扼反对称分量 f

19、unction xec,xoc=circevod(x)N=length(x);n=0:(N-1);xec=0.5*(x+x(mod(-n,N)+1); % mod取余，例如 mod(4,5)=4， mod(-4,5)= mod(5-4,5)=1xoc=0.5*(x-x(mod(-n,N)+1);序列的卷积定义： 序列的卷积函数 convu(x1,n1,x2,n2)function y,ny=convu(x1,n1,x2,n2) %convu通用卷积函数%y为卷积结果序列向量，ny是y的位置向量%x1和x2是有限长序列%n1和n2分别是x1和x2的位置向量nys=n1(1)+n2(1);nyf=

20、n1(end)+n2(end);y=conv(x1,x2);ny=nys:nyf;用重叠相加法计算序列线性卷积y=fftfilt(h,x,M) 其中，h是系统单位脉冲响应向量，x是输入序列向量，y是系统的输出序列向量（h与x的卷积结果），M是由用户选择的输入序列x的分段长度，缺省M时，默认输入序列x的分段长度M=512。使用方法如下：Lx=41;N=5;M=10; %Lx为序列x(n)长度hn=ones(1,N);hn1=hn,zeros(1,Lx-N); %产生h(n),其后补零是为了绘图好看n=0:L-1;xn=cos(pi*n/10)+cos(2*pi*n/5); %产生x(n)的Lx个

21、样值yn=ftfilt(hn,xn,M) %用重叠相加法计算卷积连续信号的卷积定义： 连续信号的卷积函数 sconv(f1,k1,f2,k2,p)functionf,k=sconv(f1,k1,f2,k2,p) %计算连续信号卷积f(t)=f1(t)*f2(t)%f：卷积f(t)对应的非零样值向量%k：f(t)对应的时间向量%f1：f1(t)的非零样值向量%k1：f1(t)对应的时间向量%f2：f2(t)的非零样值向量%k2：f2(t)对应的时间向量%p：抽样时间间隔f=conv(f1,f2); %计算f1(t)和f2(t)的卷积和f(t)f=f*p; k0=k1(1)+k2(1); %计算f

22、(t)非零样值得起点位置 k3=length(f1)+length(f2)-2;%计算卷积和f(t)的非零样值的宽度 k=k0:p:k3*p;%确定卷积和f(t)非零样值的时间向量 subplot(2,2,1) plot(k1,f1) %在子图1绘制f1(t)的时域波形图title('f1(t)') xlabel('t') ylabel('f1(t)') grid onsubplot(2,2,2) plot(k2,f2)%在子图2绘制f2(t)的时域波形图 title('f2(t)') xlabel('t') yl

23、abel('f2(t)') grid onsubplot(2,2,3) plot(k,f);%在子图3绘制卷积和f(t)的时域波形图 h=get(gca,'position') ;h(3)=2.5*h(3); set(gca,'position',h)%将第三个子图的横坐标范围扩大为原来的2.5倍 title('f(t)=f1(t)*f2(t)')xlabel('t') ylabel('f(t)')grid on序列的循环移位函数 function y=cirshftt(x,m,N)if lengt

24、h(x)>N 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); 时域循环卷积定理：设有限长序列为和，它们的N点DFT分别为和，如果，则其IDFT为两序列的循环卷积。计算两序列的循环卷积函数 function y=circonv(x1,x2,N)%循环卷积的长度N应大于等于max(length(x1),length(x2)if length(x1)>N error('N must be

25、>= 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(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&#

26、39; %用矩阵相乘的方法得到结果信噪比设纯净信号为，噪声信号为，带噪信号为，则信噪比的定义如式（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-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

27、) 要求n1<=n0<=n2function x,n=stepseq(n0,n1,n2)n=n1:n2; x=(n-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=zeros(1,length(n); y2=y1;y1(find(n>=min(n1)&(n<=ma

28、x(n1)=1)=x1; y2(find(n>=min(n2)&(n<=max(n2)=1)=x2;y=y1+y2;%4.信号乘 y(n)=x1(n)*x2(n)function y,n=sigmult(x1,n1,x2,n2)n=min(min(n1),min(n2):max(max(n1),max(n2);y1=zeros(1,length(n); y2=y1;y1(find(n>=min(n1)&(n<=max(n1)=1)=x1;y2(find(n>=min(n2)&(n<=max(n2)=1)=x2;y=y1.*y2;%5.

29、移位 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('kroneckerDelta(n,0)');%定义单位冲激单位阶跃函数 heaviside(t)使用方法如下：syms theaviside(t)录音函数function record(n,fs,file

30、) %录音函数%设置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); %保存为文件为file fprintf('The recording is saved as ');

31、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中的存储数据写入文件名为filename的文件。文件名输入是一个用单引号括起来的字符串。默认数据采样率

32、为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赫兹。单通道(单声道)音频,y是一个米长、1列向量,m是音频样本的数量。如果你的系统支

33、持立体声播放,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必须是一个两列矩阵。使用方法如下：y,fs=wavread('sound.wav

34、');%读取音频文件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

35、.wav y1=wavread('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波