1、學(xué)習(xí)好資料 歡迎下載 安裝好 MATLAB 2022 后再安裝目錄下點(diǎn)擊 會(huì)顯現(xiàn) 查找安裝程序類時(shí)出錯(cuò), 查找類時(shí)顯現(xiàn)反常 的錯(cuò)誤提示； 該錯(cuò)誤的解決方法是進(jìn)入安裝目錄下的 bin 文件夾雙擊 matlab.exe 對(duì)安裝程序進(jìn)行激活；這是可以對(duì)該 創(chuàng) 建桌面快捷方式,以后運(yùn)行程序是直接雙擊該快捷方式即可； 信號(hào)運(yùn)算 1, 信號(hào)加 MATLAB實(shí)現(xiàn) : x=x1+x2 2, 信號(hào)推遲 yn=xn-k信3, 4, 5, 號(hào)乘x=x1.*x2 信號(hào)變化幅度 y=k*x信號(hào)翻 轉(zhuǎn)y=fliplrx 信號(hào)采樣和 6, 數(shù)學(xué)描述： y= n 2 xn n n1 7, 信號(hào)采樣積 MATLAB 實(shí)現(xiàn)： y

2、=sumxn1:n2 數(shù)學(xué)描述： y n 2 xn MATLAB 實(shí)現(xiàn)： n n1 y=prodxn1:n2 8, 信號(hào)能量 數(shù)學(xué)描述： Ex n2 | xn | MATLAB實(shí)現(xiàn)： Ex=sumabsx2第 1 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 9, 信號(hào)功率 數(shù)學(xué)描述： P x 1N 1 2 | xn |Nn 0 MATLAB實(shí)現(xiàn)： Px=sumabsx2/N MATLAB 窗函數(shù) 矩形窗 w=boxcarn巴特利特窗 w=bartlettn三角窗 w=triangn布萊克曼窗 w=blackmannw=blackmann,sflag海明窗 w=haimingn其取 W=haiming

3、n,sflag sflag 用來把握窗函數(shù)首尾的兩個(gè)元素 值, 值為 symmetric, periodic漢寧窗 w=hanningn凱塞窗 w=Kaisern,beta ,beta 用于把握旁瓣的高度； n 確定時(shí), beta 越大, 其頻譜的旁瓣越小,但主瓣寬度相應(yīng)增加；當(dāng) 變化,其旁瓣高度不變； beta 確定時(shí), n 發(fā)生 切比雪夫窗： 主瓣寬度最小, 具有等波紋型, 切比雪夫窗在邊沿的采 樣點(diǎn)有尖峰； W=chebwinn,r第 2 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 數(shù)字濾波器的特性分析 1,脈沖響應(yīng)： impz 函數(shù) 調(diào)用方式： （ 1）h,t=impzb,a:返回參數(shù) h

4、是脈沖響應(yīng)的數(shù)據(jù), t 是脈沖響應(yīng)的 時(shí)間間隔； （ 2）h,t=impzb,a,N: N 用來指定脈沖信號(hào)的長(zhǎng)度； （ 3）h,t=impzb,a,n,Fs:Fs 用來指定脈沖信號(hào)的采樣頻率 4 h,t=impzb,a,Fs:不再指定指定脈沖信號(hào)的長(zhǎng)度； 例： b,a=butter4,0.05; impzb,a,100; 2,頻率響應(yīng)（幅頻響應(yīng)和相頻響應(yīng)） （ 1）數(shù)字濾波器頻率響應(yīng)： freqz 函數(shù) 調(diào)用方式： h,w=freqzb,a,n: 返回濾波器的 n 點(diǎn)復(fù)頻率響應(yīng), b,a 分別是濾波器系數(shù)的分子和 分母向量； h 是復(fù)頻率響應(yīng), w 是頻率點(diǎn)； h,w=freqzb,a,n,

5、 whole:接受單位圓上的 n 個(gè)點(diǎn)； 第 3 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 h= freqz b,a,w h,f=freqzb,a,n,fs h= freqz b,a,f,fs 2模擬濾波器頻率響應(yīng)： freqs 函數(shù) 調(diào)用方式： h=freqsb,a,w：w 指定頻率點(diǎn)的復(fù)頻率響應(yīng) h,w=freqsb,a,n:用 n 指定進(jìn)行復(fù)頻率響應(yīng)的采樣點(diǎn)數(shù) 例： 1; 1; w=logspace-1,1;freqsb,a,w;10 e ud it n g a M 10 -1 10 0 10 1 10 Frequency rad/s 50 s ere g ed e s a hP 00 10

6、 1 10 -50 -100 -1 10 Frequency rad/s 3,幅頻和相頻 y=absx:運(yùn)算 x 各元素的確定值；當(dāng) x 為一個(gè)復(fù)數(shù)時(shí),運(yùn)算 x第 4 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 的復(fù)數(shù)模； Y=anglex: 運(yùn)算 x 向量各元素的復(fù)數(shù)相位值,單位為弧度； 功率譜估（） 一, 隨機(jī)信號(hào)處理基礎(chǔ) 1, mean 函數(shù) 調(diào)用方式： （ 1） y=meanX：當(dāng) X 為向量時(shí),此函數(shù)結(jié)果為 X 的均值；當(dāng) X 為矩陣時(shí), 函數(shù)結(jié)果為一個(gè)行向量,其元素分別為矩陣每列元素的均值； （ 2） （ 3） 2, y=meanX,dim: dim=1 時(shí),函數(shù)結(jié)果為一個(gè) 行向量 ,其

7、元素分別為矩陣 每列元素的均值 ； din=2 時(shí),函數(shù)結(jié)果為一個(gè) 列向量 ,其元素分別為矩陣 每行元素的均值 ； 協(xié)方差： cov 函數(shù) 調(diào)用方式 （ 1） y=covX:當(dāng) X 為向量 時(shí),函數(shù)返回結(jié)果為 X 的方差 ；當(dāng) X 為矩陣 時(shí), 就它的每一列相當(dāng)于一個(gè)變量,函數(shù)返回結(jié)果為該矩陣的 列與列之間的 協(xié)方差矩陣 ,diagcovX是該矩陣每一個(gè)列向量的方差； （ 2） y=covX,Y:相當(dāng)于 covX:y: ,運(yùn)算兩個(gè)等長(zhǎng)度向量的互協(xié)方差矩陣； 例如： X=1 2 3 4;5 6 7 8;2 4 7 8;1 0 2 3;4 5 6 7； 3, A=covX； % 運(yùn)算協(xié)方差 % 運(yùn)

8、算互協(xié)方差 B=covX1,:,X3,:;相關(guān)函數(shù)估量 （ 1） xcorr 函數(shù)：相關(guān)函數(shù)估量 C=xcorrA,B:當(dāng)和為長(zhǎng)度為（）的向量時(shí),返回結(jié)果為長(zhǎng) 度為 2M-1 的相互關(guān)函數(shù)序列；當(dāng) A 和 B 長(zhǎng)度不同時(shí),就要對(duì)長(zhǎng)度小的進(jìn)行補(bǔ) 零；假如 A 為列向量,就 C 也為列向量,假如 A 為行向量,就 C 也為行向量； C=xcorrA: C=xcorrA:估量向量 A 的自相關(guān)函數(shù)； 當(dāng) A 為 M*N 的矩陣時(shí),返回結(jié)果為（ 2M-1 ）行,N2 列的 矩陣,該矩陣的列是由矩陣 A 全部列之間的相互關(guān)函數(shù)構(gòu)成； C=xcorr, scaleopt: 參數(shù) scaleopt 用來指定

9、相關(guān)函數(shù)估量所接受 估量方式,即 biased:有偏估量方式 unbiased: 無(wú)偏估量方式 coeff: none： 對(duì)序列進(jìn)行歸一化處理 運(yùn)算序列的非歸一化相關(guān) （ 2） xcov 函數(shù)：協(xié)方差函數(shù)估 量 （ 3） 相關(guān)系數(shù)估量運(yùn)算 第 5 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 Corrcoef 函數(shù)：運(yùn)算序列的相關(guān)系數(shù) 二, 經(jīng)典功率譜估量方法 1, 直接法（周期圖法） Periodogram 函數(shù)：功率譜估量 Pxx=periodogramx:返回向量 x 的功率譜估量向量 Pxx.Pxx=periodogramx,window: 參數(shù) window 用來指定所接受的窗函數(shù) Pxx,

10、w=periodogramx,window,NFFT: 如 x 為實(shí)信號(hào), NFFT 為偶數(shù),就 Pxx 的長(zhǎng)度為（ NFFT/2+1 ）； 如 x 為實(shí)信號(hào), NFFT 為奇數(shù),就 Pxx 的長(zhǎng)度為（ NFFT+1 ） /2； 如 x 為復(fù)信號(hào),就 如 x 為實(shí)信號(hào),就 如 x 為復(fù)信號(hào),就 Pxx 的長(zhǎng)度為 NFFT ； w 的范疇為 0,Pi; w 的范疇為 0,2*Pi;Pxx,f=periodogramx,window,NFFT,Fs: 同時(shí)返回和估量 PSD 的位置一一對(duì)應(yīng)的線性頻率 f, 參數(shù) Fs 為采樣頻率； 如 x 為實(shí)信號(hào),就 如 x 為復(fù)信號(hào),就 f 的范疇為 0,Fs

11、/2; f 的范疇為 0,Fs;例 1：接受 periodogram 函數(shù)來運(yùn)算功率譜 Fs=2022; NFFT=1024; n=0:1/Fs:1; x=sin2*pi*100*n+4*sin2*pi*500*n+randnsizen; window=boxcarlengthx; periodogramx,window,NFFT,Fs;第 6 頁(yè),共 18 頁(yè)10 學(xué)習(xí)好資料 歡迎下載 Periodogram Power Spectral Density Estimate 0-10 z /HB d y c nue q e fr r/ e w oP -20 -30 -40 -50 -60 -7

12、0 -80 00.1 Frequency kHz 1例 2,利用 FFT 直接法運(yùn)算上面噪聲信號(hào)的功率譜 Fs=2022; nFFT=1024; n=0:1/Fs:1; x=sin2*pi*100*n+4*sin2*pi*500*n+randnsizen; X=fftx,nFFT; Pxx=absX.2/lengthn; t=0:roundnFFT/2-1; k=t*Fs/nFFT; p=10*log10Pxxt+1; plotk,p;40 30 20 10 0-10 -20 -30 2, -40 0100 200 300 400 500 600 700 800 900 1000 間接法 第

13、7 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 40 35 30 25 20 15 10 500100 200 300 400 500 600 700 800 900 1000 Fs=2022; nFFT=1024; n=0:1/Fs:1; x=sin2*pi*100*n+4*sin2*pi*500*n+randnsizen; Cx=xcorrx,unbiased; Cxk=fftCx,nFFT; Pxx=absCxk; t=0:roundnFFT/2-1; k=t*Fs/nFFT; p=10*log10Pxxt+1;plotk,p;3, 協(xié)方差法估量 自回來功率譜估量的協(xié)方差方法 Pcov 函自回

14、來功率譜估量的改進(jìn)的協(xié)方差方法 數(shù) Pmcov 函 數(shù) 應(yīng)用實(shí)例：比較兩種方法在噪聲信號(hào)的功率譜估量中的成效,發(fā)覺兩種方法基 本相同； Fs=1000; h=fir120,0.3; r=randn1024,1; x=filterh,1,r; p1,f=pcovx,20,Fs; p2,f=pmcovx,20,Fs; Pxx1=10*log10p1; Pxx2=10*log10p2; plotf,Pxx1,r:,f,Pxx2,g-; ylabel幅值 :dB; xlabel功率譜估量 ;legend協(xié)方差方法 , 改進(jìn)的協(xié)方差方法 ;第 8 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 -20 協(xié) 方 差

15、 方 法 -30 改 進(jìn) 的 協(xié) 方 差 方 法 -40 -50 -60 -70 -80 -90 -100 050 100 150 200 250 300 350 400 450 500 三, 現(xiàn)在譜估量的非參數(shù)方法 1, MTM （ Multitaper ）法估量 時(shí)間 -帶寬的乘積 NW ,窗口數(shù) =2*NW-1 Pmtm 函數(shù)：實(shí)現(xiàn) Multitaper 法的功率譜估量 調(diào)用方式 （ 1） Pxx=pmtmx: 用 Multitaper 法對(duì)離散時(shí)間信號(hào) x 進(jìn)行功率譜估量, 如 x 為實(shí)數(shù),就返回結(jié)果為“單邊” 功率譜,如 x 為復(fù)數(shù),就返回結(jié)果為 “雙邊” 功率譜； （ 2） Pxx

16、=pmtmx,NW:（ 3） Pxx=pmtmx,NW,NFFT: 參數(shù) NFFT 用來指定 FFT 運(yùn)算所接受的點(diǎn)數(shù)； 如 x 為實(shí)信號(hào), NFFT 為偶數(shù),就 Pxx 的長(zhǎng)度為（ NFFT/2+1 ）； 如 x 為實(shí)信號(hào), NFFT 為奇數(shù),就 Pxx 的長(zhǎng)度為（ NFFT+1 ） /2； 如 x 為復(fù)信號(hào),就 Pxx 的長(zhǎng)度為 NFFT ； NFFT 默認(rèn)值為 256（ 4） Pxx,w=pmtm：輸出參數(shù) w 和估量 PSD 的位置一一對(duì)應(yīng)的歸一化 角頻率,單位 rad/sample,范疇如下： 如 x 為實(shí)信號(hào),就 如 x 為復(fù)信號(hào),就 w 的范疇為 0,Pi; w 的范疇為 0,2

17、*Pi;第 9 頁(yè),共 18 頁(yè)P(yáng)xx,f=pmtm 同時(shí)返回和估量 學(xué)習(xí)好資料 歡迎下載 , Fs: PSD 的位置一一對(duì)應(yīng)的線性頻率 f, 參數(shù) Fs 為采樣頻率； 如 x 為實(shí)信號(hào),就 f 的范疇為 0,Fs/2; 如 x 為復(fù)信號(hào),就 f 的范疇為 0,Fs;（ 5） Pxx,f=pmtm, Fs,mehod:Method 有： adapt:Thomson 自適應(yīng)非線性組合算法,默認(rèn) 值 Unity: 相同加權(quán)的線性組合 （ 6） Pxx,Pxxc,f=pmtm Eigen:特點(diǎn)值加權(quán)的線性組合 , Fs,mehod,p: P:01,指定 PSD 的置信度； Pxxc:,1是置信度區(qū)間

18、的下部分的數(shù)值, 是置信度區(qū)間的上部分的數(shù)值； Pxxc:,2（ 7） Pxx,Pxxc,f=pmtmx, dpss_params,NFFT,Fs,mehod,p:（ 8） =pmtm, twoside/oneside：應(yīng)用說明：利用 Multitaper 進(jìn)行 PSD 估量,并比較 NW取不同數(shù)值時(shí)的結(jié)果； Fs=1000; nFFT=1024; t=0:1/Fs:1; x=sin2*pi*200*t+randnsizet; P1,f=pmtmx,2,nFFT,Fs; P2,f=pmtmx,4,nFFT,Fs;P3,f=pmtmx,10,nFFT,Fs; Pxx1=10*log10P1; P

19、xx2=10*log10P2; Pxx3=10*log10P3; subplot3,1,1; plotf,Pxx1; xlabelNW=2; subplot3,1,2; plotf,Pxx2; xlabelNW=4; subplot3,1,3; plotf,Pxx3; xlabelNW=10;第 10 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 0-20 -40 050 100 150 200 250 300 350 400 450 500 NW=2 0-20 -40 050 100 150 200 250 300 350 400 450 500 NW=4 0-20 -40 050 100 150

20、200 250 300 350 400 450 500 NW=10 可以看出： NW數(shù)值越大, 曲線越平滑, 說明方差比較??； 但同時(shí)看到波峰變寬, 這說明頻譜泄露增大了； 2,MUSIC 法估量 Pmusic 函數(shù)：實(shí)現(xiàn) MUSIC 法的功率譜運(yùn)算調(diào)用方式 1 s=pmusicx,p: 用 MUSIC 法對(duì)離散時(shí)間信號(hào) x 進(jìn)行功率譜估量； p 是信號(hào) x中包含的復(fù)數(shù)正弦波信號(hào)的個(gè)數(shù),假如 x 是一個(gè)數(shù)據(jù)矩陣,就對(duì)矩陣的每一 列都進(jìn)行功率譜估量；留意：為了返回實(shí)信號(hào)的全部功率譜值,需要使用另 一個(gè)輸入?yún)?shù) whole.（ 3） s=pmusicr,p,corr: r來指定自相關(guān)矩陣； 對(duì)于實(shí)

21、信號(hào)而言, 默認(rèn)情形下 Pmusic 返回功率譜估量的一半； 對(duì)于復(fù)信號(hào), 返回全部的功率譜估量（ 4） 值； （ 5） s=pmusicx,p,NFFT: Pxx,w=pmusic （ 6） Pxx,w=pmusic ,Fs第 11 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 （ 7） s,f=pmusic ,NW,noverlapNW 默認(rèn)值為 2*p, 參數(shù) noverlap 的 默認(rèn)值為 NW-1（ 8） s,w,v,e=pmusic ： 輸出參數(shù) v 為以矩陣,其列是與噪聲子空間一 一對(duì)應(yīng)的特點(diǎn)值所組成的向量；而 e 為相關(guān)矩陣的特點(diǎn)值向量； 應(yīng)用說明 例：用 MUSIC 法進(jìn)行 PSD

22、估量,結(jié)果如以下圖： randnstate; n=0:99;s=expi*pi/2*n+2*expi*pi/4*n+expi*pi/3*n+randn1,100;X=corrmtxs,7,mod; % 使用改進(jìn)的協(xié)方差方法估量相互關(guān)矩陣 s,w,e,v=pmusicX,4,whole;pmusicX,4,whole;70 Pseudospectrum Estimate v ia MUSIC 60 50 40 dB r e w o P 30 20 10 0-10 -20 0Normalized Frequency 1rad/sample 3,特點(diǎn)向量（ AV ）法估量（也是基于矩陣特點(diǎn)分解的一種

23、功率譜估量的非參數(shù) 方法,它主要適用于混有在噪聲的正弦信號(hào)的功率譜估量） Peig 函數(shù)：實(shí)現(xiàn)特點(diǎn)向量法的功率譜估量 調(diào)用方式 （ 1）s=peigx,p: 用特點(diǎn)向量法對(duì)離散時(shí)間信號(hào) x 進(jìn)行功率譜估量； p 是信號(hào) x 中包含的復(fù)數(shù)正弦波信號(hào)的個(gè)數(shù)； 2s=peigr,p,corr3s=peigx,p,nFFT4Pxx,w=peig 5 Pxx,f=peig 6 ,Fs, NW ,noverlaps,f=peig第 12 頁(yè),共 18 頁(yè)7s,w,v,e=peig學(xué)習(xí)好資料 歡迎下載 8peig應(yīng)用說明 例：用特點(diǎn)向量法進(jìn)行 PSD 估量,結(jié)果如下圖所示： randnstate,1; n=

24、0:99; s=expi*pi/2*n+2*expi*pi/4*n+expi*pi/3*n+randn1,100; X=corrmtxs,7,mod; s,w,e,v=pmusicX,4,whole; peigX,4,whole;70 Pseudospectrum Estimate via Eigenvector Method 60 50 40 B d r e w o P 30 20 10 0-10 -20 -30 01rad/sample Normalized Frequency MUSIC 算法 m=sqrt-1; delta=0.101043; a1=-0.850848; sample=

25、32; %number of sample spot p=10; %number of sample spot in coef method; 第 13 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 f1=0.05; f2=0.40; f3=0.3; fstep=0.01; fstart=-0.5;fend=0.5; f=fstart:fstep:fend; nfft=fend-fstart/fstep+1; %un=urn+juin urn= normrnd0,delta/2,1,sample;uin= normrnd0,delta/2,1,sample;un=urn+m*uin; %. zn fo

26、r n=1:sample-1 zn1=un1; znn+1=-a1*znn+unn+1; end %. xn for n=1:sample xnn=2*cos2*pi*f1*n-1+2*cos2*pi*f2*n-1+2*cos2*pi*f3*n-1+ sqrt2*realznn; end %* x=xn; for k=0:1:sample-1 s=0; for n=1:sample-k, s=s+conjxn,1*xn+k,1; end rxx1,k+1=1/sample*s; end Rx=zerossample,sample; Rx=toeplitzrxx1,1:32; U,S,V=svd

27、Rx; Pmusicf=zeros1,1/fstep+1; ei=zeros1,sample; for i=1:lengthf %calculate the value of rxx for j=1:sample ei1,j=exp-2*pi*j-1*fi*m; end ; sum=0; for k=7:sample sum=sum+absei*V:,k2; 第 14 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 end Pmusicf1,i=10*log101/sum;end figure plotf,Pmusicf; 10 50-5 -10 -15 -20 0第 15 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料

28、歡迎下載 0-10 B d e d u it n g am -20 -30 -40 -50 -60 -90 -60 -30 030 60 90 angle degree derad = pi/180; % deg - rad radeg = 180/pi; twpi = 2*pi; kelm = 8;%陣列數(shù)量 dd = 0.5;% spaced=0:dd:kelm-1*dd; %iwave = 3;% number of DOAtheta = 10 30 60;角度 snr = 10;% input SNR dBn = 500;%A=exp-j*twpi*d.*sintheta*derad;

29、% direction matrix S=randniwave,n; X=A*S; X1=awgnX,snr,measured; Rxx=X1*X1/n; InvS=invRxx; % EV,D=eigRxx;% EVA=diagD; EVA,I=sortEV A; EVA=fliplrEV A; EV=fliplrEV:,I;% MUSIC for iang = 1:361第 16 頁(yè),共 18 頁(yè)學(xué)習(xí)好資料 歡迎下載 angleiang=iang-181/2; phim=derad*angleiang; a=exp-j*twpi*d*sinphim.; L=iwave; En=EV:,L+1:kelm; SPiang=a*a/a*En*En*a;end% SP=absSP; SPmax=maxSP; SP=10*log10SP/SPmax; h=plotangle,SP; seth,Linewidth,2 xla