matlab信號處理相關(guān)函數(shù)程序

1、序列的傅里葉變換及其逆變換定義： 其幅度特性為，在Matlab中采用abs函數(shù)；相位特性為，在Matlab中采用angle函數(shù)。為了方便，考慮在兩個周期，例如中2M+1個均勻頻率點上計算FT，并且觀察其周期性和對稱性。為此給出function文件如下，求解FT變換：functionX,w=ft(x,n,k)%X：序列x(n)的傅里葉變換%w：X的自變量%x：要進行傅里葉變換的序列x(n)%n：序列x(n)的位置向量%k：求和區(qū)間w=(pi/abs(max(k)/2)*k;X=x*(exp(-1i*pi/abs(max(k)/2).(n*k);使用方法如下：n=-5:5;%序列區(qū)間x=(-0.9

2、).n;%序列表達式k=-200:200;%求和區(qū)間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定義： 離散傅里葉變換的的性質(zhì)：（1）DFT的共軛對稱性。若，則：， 。（2）實序列DFT的性質(zhì)。若為實序列，則其離散傅里葉變換為共軛對稱，即。（3）實偶序列DFT的性質(zhì)。若為實偶序列，則其離散傅里葉變換為實偶對稱，即。（4）實

5、奇序列DFT的性質(zhì)。若為實奇序列，則其離散傅里葉變換為純虛奇對稱，即。離散傅立葉變換函數(shù) 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);離散傅立葉反變換函數(shù) 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); %產(chǎn)生矩形序列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、序列產(chǎn)生及其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、換函數(shù) fftXk=fft(xn,N) xn為時域序列向量，N是DFT變換區(qū)間長度。當(dāng)N大于xn的長度時，fft函數(shù)自動在xn后面補零。函數(shù)返回xn的N點DFT變換結(jié)果向量Xk。當(dāng)N小于xn的長度時，fft函數(shù)計算xn的前面N個元素構(gòu)成的N長序列的N點DFT，忽略xn后面的元素。當(dāng)x為矩陣時，Xk為矩陣xn的每一列的FFT ,fft 函數(shù)按類似的方法處理列長度。當(dāng)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;%產(chǎn)生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)連續(xù)信號的傅里葉變換及其逆變換定義： Fw=fourier(ft,t,w)%求時域函數(shù)ft的Fourier變換FwFt=ifourier(Fw,w,t)%求頻域函數(shù)Fw的Fourie

10、r逆變換ft使用方法如下：syms t wut=heaviside(t);%單位階躍函數(shù)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逆變換連續(xù)信號的單邊拉普拉斯變換及其逆變換定義：Fs=laplace(ft,t,s)%求時域函數(shù)ft的Laplace變換Fsft=ilaplace(Fs,s,t)%求頻域函數(shù)Fs的Laplace逆變換ft使用方法如下：syms t ssyms a positiveft=dirac(t);%單位沖激函數(shù)Fs=laplace(ft,t,s);%拉普拉斯變換Inv_ft=ifourier(Fs,s,t);%拉普拉斯逆變換線性常系數(shù)差分方程的遞推求解 求差分方程的零狀態(tài)響應(yīng)和全響應(yīng)函數(shù)yn=filter(B,A,xn) 計算系統(tǒng)對輸入信號向量xn的

12、零狀態(tài)響應(yīng)輸出信號向量yn，yn與xn長度相等，其中，B和A是差分方程的系數(shù)向量，即 ,其中，如果，則filter用 對系數(shù)向量B和A歸一化。yn=filter(B,A,xn,xi) 計算系統(tǒng)對輸入信號向量xn的全響應(yīng)輸出信號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，調(diào)用時可缺省xs。使用方法如下：%調(diào)用filter解差分方程y(n)-ay(n-1)=x(

13、n)a=0.8;ys=1; %設(shè)差分方程系數(shù)a=0.8，初始狀態(tài)：y(-1)=1xn=1,zeros(1,30); %x(n)=單位脈沖序列，長度N=31B=1;A=1,-a; %差分方程系數(shù)xi=filtic(B,A,ys); %由初始條件計算等效初始條件的輸入序列xiyn=filter(B,A,xn,xi); %調(diào)用filter解差分方程，求系統(tǒng)輸出函數(shù)y(n)n=0:length(yn)-1;stem(n,yn,.)xlabel(n);ylabel(y(n)離散系統(tǒng)函數(shù)定義： 系統(tǒng)的因果性：其單位脈沖響應(yīng)h(n)是因果序列，滿足n0時，h(n)=0或 其系統(tǒng)函數(shù)H(z)的收斂域包含無窮點

14、系統(tǒng)的穩(wěn)定性：其單位脈沖響應(yīng)h(n)絕對可和，即 或 其系統(tǒng)函數(shù)H(z)的收斂域包含單位圓系統(tǒng)穩(wěn)定性判定函數(shù) stab(A)function stab(A)%stab：系統(tǒng)穩(wěn)定性判定函數(shù)，A是H(z)的分母多項式系數(shù)向量disp(系統(tǒng)極點為：)P=roots(A)%求H(z)的極點，并顯示disp(系統(tǒng)極點模的最大值為：)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、; 時域循環(huán)卷積定理：設(shè)有限長序列為和，它們的N點DFT分別為和，如果，則其IDFT為兩序列的循環(huán)卷積。計算兩序列的循環(huán)卷積函數(shù) function y=circonv(x1,x2,N)%循環(huán)卷積的長度N應(yīng)大于等于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); %該語句的功能相當(dāng)于序列翻褶，延拓，取主值序列H=zeros(N,N);for n=1:1:N %得到x3序列的循環(huán)移位矩陣 H(n,:)=cirshftt(x3,n-1,N); endy=x1*H %用矩陣相乘的方法得到結(jié)果信噪比設(shè)純凈信號為，噪聲信號為，帶噪信號為，則信噪比的定義如式（4-2）所示，單位為dB。 （4-2）function y=snr(x0,xn)%計算信噪比，單位dB%x0是純凈信號,xn是帶噪信號p0=sum(abs(x0).2); %sum函數(shù)為求和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函數(shù)：找出非零元素的索引號%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);特殊的連續(xù)函數(shù)單位沖激函數(shù) dirac(t)使用方法如下：syms t ndirac(t)delta=sym(kronec

19、kerDelta(n,0);%定義單位沖激單位階躍函數(shù) heaviside(t)使用方法如下：syms theaviside(t)錄音函數(shù)function record(n,fs,file) %錄音函數(shù)%設(shè)置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的抽樣數(shù)據(jù)。如果不包括一個擴展文件名,wavread默認文件后綴為. wav。y, Fs = wavread(filename) 返回的采樣率(Fs)赫茲用于編碼中的數(shù)據(jù)文件。y, Fs, nbits = wavread(filename) 返回每個抽樣的比特數(shù)(nbits)。寫入音頻wavwrite(y,filename) 把在變量y中的存儲數(shù)據(jù)寫入文件名為filen

21、ame的文件。文件名輸入是一個用單引號括起來的字符串。默認數(shù)據(jù)采樣率為8000 Hz和被認為是16位。每一列代表一個單獨的數(shù)據(jù)通道。因此,立體數(shù)據(jù)應(yīng)該指定為一個矩陣有兩個列。wavwrite(y,Fs,filename) 把在變量y中的存儲數(shù)據(jù)寫入文件名為filename的文件。數(shù)據(jù)的采樣率為Fs赫茲和被認為是16位。wavwrite(y,Fs,N,filename) 把在變量y中的存儲數(shù)據(jù)寫入文件名為filename的文件。數(shù)據(jù)采樣率為Fs赫茲和N位,其中N是8,16、24、32。播放音頻sound(y,Fs) 以采樣率Fs播放音頻信號y。如果你不指定采樣率,默認為8192赫茲。單通道(單聲

22、道)音頻,y是一個米長、1列向量,m是音頻樣本的數(shù)量。如果你的系統(tǒng)支持立體聲播放,y是一個m-by-2矩陣,第一列對應(yīng)于左聲道,和第二列對應(yīng)于右聲道。聲音功能假設(shè)y包含-1和1之間的浮點數(shù),和剪輯范圍之外的值。sound(y,Fs,bits) 指定的位深度(即精度)的樣本值。位深度的可能的值取決于可用的音頻硬件系統(tǒng)上。大多數(shù)平臺支持8位或16位的深度。如果你不指定位,聲音功能在一個8位深度。wavplay(y,Fs) 播放存儲的音頻信號向量y。Fs是整數(shù)采樣率，單位是Hz。Fs的默認值是11025赫茲。wavplay只支持1 或2聲道音頻(單聲道或立體聲)信號。在立體聲,y必須是一個兩列矩陣。

23、使用方法如下：y,fs=wavread(sound.wav);%讀取音頻文件L=length(y);%采樣點數(shù)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個點進行計