2022年matlab与信号处理知识点 2.pdf

《2022年matlab与信号处理知识点 2.pdf》由会员分享，可在线阅读，更多相关《2022年matlab与信号处理知识点 2.pdf（18页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、名师精编优秀资料安装好 MATLAB 2012 后再安装目录下点击setup.exe 会出现查找安装程序类时出错， 查找类时出现异常 的错误提示。该错误的解决方法是进入安装目录下的bin 文件夹双击 matlab.exe 对安装程序进行激活。这是可以对该matlab.exe 创建桌面快捷方式，以后运行程序是直接双击该快捷方式即可。信号运算1、信号加MATLAB 实现: x=x1+x2 2、信号延迟y(n)=x(n-k) 3、信号乘x=x1.*x2 4、信号变化幅度y=k*x 5、信号翻转y=fliplr(x) 6、信号采样和数学描述： y=21)(nnnnx MATLAB 实现： y=sum(

2、x(n1:n2) 7、信号采样积数学描述：21)(nnnnxy MATLAB实现： y=prod(x(n1:n2) 8、信号能量数学描述：nxnxE2| )(|MATLAB 实现：Ex=sum(abs(x)2) 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 1 页，共 18 页 - - - - - - - - - 名师精编优秀资料9、信号功率数学描述：102|)(|1PNnxnxNMATLAB 实现： Px=sum(abs(x)2)/N MATLAB窗函数矩形窗w

3、=boxcar(n) 巴特利特窗w=bartlett(n) 三角窗w=triang(n) 布莱克曼窗w=blackman(n) w=blackman(n,sflag) 海明窗w=haiming(n) W=haiming(n,sflag) sflag 用来控制窗函数首尾的两个元素值，其取值为 symmetric 、periodic 汉宁窗w=hanning(n) 凯塞窗w=Kaiser(n,beta) ,beta 用于控制旁瓣的高度。 n 一定时， beta 越大，其频谱的旁瓣越小，但主瓣宽度相应增加；当beta 一定时， n 发生变化，其旁瓣高度不变。切比雪夫窗： 主瓣宽度最小， 具有等波纹型

4、， 切比雪夫窗在边沿的采样点有尖峰。W=chebwin(n,r) 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 2 页，共 18 页 - - - - - - - - - 名师精编优秀资料数字滤波器的特性分析1、脉冲响应： impz 函数调用方式：（1）h,t=impz(b,a): 返回参数 h 是脉冲响应的数据， t 是脉冲响应的时间间隔。（2）h,t=impz(b,a,N): N用来指定脉冲信号的长度。（3）h,t=impz(b,a,n,Fs):Fs用来指定脉冲

5、信号的采样频率(4) h,t=impz(b,a,Fs): 不再指定指定脉冲信号的长度。例：b,a=butter(4,0.05); impz(b,a,100); 2、频率响应（幅频响应和相频响应）（1）数字滤波器频率响应：freqz 函数调用方式：h,w=freqz(b,a,n): 返回滤波器的n 点复频率响应， b,a 分别是滤波器系数的分子和分母向量。 h 是复频率响应， w 是频率点。h,w=freqz(b,a,n, whole ):采用单位圆上的n 个点。名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - -

6、- - - - - - - - - - 第 3 页，共 18 页 - - - - - - - - - 名师精编优秀资料h= freqz (b,a,w) h,f=freqz(b,a,n,fs) h= freqz (b,a,f,fs) (2)模拟滤波器频率响应：freqs 函数调用方式：h=freqs(b,a,w) ：w 指定频率点的复频率响应h,w=freqs(b,a,n): 用 n 指定进行复频率响应的采样点数例： b=0.3 0.6 1; a=1 0.5 1; w=logspace(-1,1);freqs(b,a,w); 10-1100101-100-50050Frequency (rad/

7、s)Phase(degrees)10-110010110-0.5100.2Frequency (rad/s)Magnitude3、幅频和相频y=abs(x): 计算 x 各元素的绝对值。当x 为一个复数时，计算x名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 4 页，共 18 页 - - - - - - - - - 名师精编优秀资料的复数模。Y=angle(x): 计算 x 向量各元素的复数相位值，单位为弧度。功率谱估（）一、随机信号处理基础1、mean 函数调用方

8、式：（1） y=mean(X) ：当 X 为向量时，此函数结果为X 的均值；当 X 为矩阵时，函数结果为一个行向量，其元素分别为矩阵每列元素的均值。（2） y=mean(X,dim): （3） dim=1 时，函数结果为一个 行向量 ，其元素分别为矩阵 每列元素的均值 。din=2 时，函数结果为一个 列向量，其元素分别为矩阵 每行元素的均值 。2、协方差： cov 函数调用方式（1） y=cov(X): 当 X 为向量时，函数返回结果为X 的方差 ；当 X 为矩阵时，则它的每一列相当于一个变量，函数返回结果为该矩阵的列与列之间的协方差矩阵 ，diag(cov(X) 是该矩阵每一个列向量的方差

9、。（2） y=cov(X,Y): 相当于 cov(X(:)y(:) ，计算两个等长度向量的互协方差矩阵。例如： X=1 2 3 4;5 6 7 8;2 4 7 8;1 0 2 3;4 5 6 7；A=cov(X) ；%计算协方差B=cov(X(1,:),X(3,:); % 计算互协方差3、相关函数估计（1） xcorr 函数：相关函数估计C=xcorr(A,B): 当和为长度为（）的向量时，返回结果为长度为 2M-1 的互相关函数序列；当A 和 B 长度不同时，则要对长度小的进行补零；如果 A 为列向量，则 C 也为列向量，如果A 为行向量，则 C 也为行向量。C=xcorr(A): 估计向量

10、 A 的自相关函数。C=xcorr(A): 当 A 为 M*N 的矩阵时，返回结果为（ 2M-1）行、N2 列的矩阵，该矩阵的列是由矩阵A 所有列之间的互相关函数构成。C=xcorr( ， scaleopt): 参数 scaleopt 用来指定相关函数估计所采用的估计方式，即biased: 有偏估计方式unbiased: 无偏估计方式coeff: 对序列进行归一化处理none：计算序列的非归一化相关（2） xcov 函数：协方差函数估计（3） 相关系数估计计算名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - -

11、 - - - - - - - - - 第 5 页，共 18 页 - - - - - - - - - 名师精编优秀资料Corrcoef 函数：计算序列的相关系数二、 经典功率谱估计方法1、直接法（周期图法）Periodogram 函数：功率谱估计Pxx=periodogram(x): 返回向量 x 的功率谱估计向量Pxx. Pxx=periodogram(x,window): 参数 window 用来指定所采用的窗函数Pxx,w=periodogram(x,window,NFFT): 若 x 为实信号， NFFT 为偶数，则 Pxx 的长度为（ NFFT/2+1） ；若 x 为实信号， NFFT

12、 为奇数，则 Pxx 的长度为（ NFFT+1）/2；若 x 为复信号，则 Pxx 的长度为 NFFT ；若 x 为实信号，则 w 的范围为 0,Pi; 若 x 为复信号，则 w 的范围为 0,2*Pi; Pxx,f=periodogram(x,window,NFFT,Fs): 同时返回和估计PSD 的位置一一对应的线性频率f, 参数 Fs 为采样频率。若 x 为实信号，则 f 的范围为 0,Fs/2; 若 x 为复信号，则 f 的范围为 0,Fs; 例 1：采用 periodogram 函数来计算功率谱Fs=2000; NFFT=1024; n=0:1/Fs:1; x=sin(2*pi*10

13、0*n)+4*sin(2*pi*500*n)+randn(size(n); window=boxcar(length(x); periodogram(x,window,NFFT,Fs); 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 6 页，共 18 页 - - - - - - - - - 名师精编优秀资料00.10.20.30.40.50.60.70.80.91-80-70-60-50-40-30-20-10010Frequency (kHz)Power/fre

14、quency(dB/Hz)Periodogram Power Spectral Density Estimate例 2、利用 FFT 直接法计算上面噪声信号的功率谱Fs=2000; nFFT=1024; n=0:1/Fs:1; x=sin(2*pi*100*n)+4*sin(2*pi*500*n)+randn(size(n); X=fft(x,nFFT); Pxx=abs(X).2/length(n); t=0:round(nFFT/2-1); k=t*Fs/nFFT; p=10*log10(Pxx(t+1); plot(k,p); 010020030040050060070080090010

15、00-40-30-20-100102030402、间接法名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 7 页，共 18 页 - - - - - - - - - 名师精编优秀资料010020030040050060070080090010000510152025303540Fs=2000; nFFT=1024; n=0:1/Fs:1; x=sin(2*pi*100*n)+4*sin(2*pi*500*n)+randn(size(n); Cx=xcorr(x,unbi

16、ased); Cxk=fft(Cx,nFFT); Pxx=abs(Cxk); t=0:round(nFFT/2-1); k=t*Fs/nFFT; p=10*log10(Pxx(t+1); plot(k,p); 3、协方差法估计自回归功率谱估计的协方差方法Pcov 函数自回归功率谱估计的改进的协方差方法Pmcov 函数应用实例：比较两种方法在噪声信号的功率谱估计中的效果，发现两种方法基本相同。Fs=1000; h=fir1(20,0.3); r=randn(1024,1); x=filter(h,1,r); p1,f=pcov(x,20,Fs); p2,f=pmcov(x,20,Fs); Pxx

17、1=10*log10(p1); Pxx2=10*log10(p2); plot(f,Pxx1,r:,f,Pxx2,g-); ylabel( 幅值:dB); xlabel( 功率谱估计 ); legend(协方差方法 , 改进的协方差方法 ); 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 8 页，共 18 页 - - - - - - - - - 名师精编优秀资料050100150200250300350400450500-100-90-80-70-60-50-40

18、-30-20协 方 差 方 法改 进 的 协 方 差 方 法三、 现在谱估计的非参数方法1、MTM （Multitaper ）法估计时间-带宽的乘积 NW ，窗口数 =2*NW-1 Pmtm 函数：实现 Multitaper法的功率谱估计调用方式（1） Pxx=pmtm(x): 用 Multitaper法对离散时间信号x 进行功率谱估计， 若x 为实数，则返回结果为“单边”功率谱，若 x 为复数，则返回结果为“双边”功率谱。（2） Pxx=pmtm(x,NW): （3） Pxx=pmtm(x,NW,NFFT): 参数 NFFT 用来指定 FFT 运算所采用的点数。若 x 为实信号， NFFT

19、为偶数，则 Pxx 的长度为（ NFFT/2+1） ；若 x 为实信号， NFFT 为奇数，则 Pxx 的长度为（ NFFT+1）/2；若 x 为复信号，则 Pxx 的长度为 NFFT ；NFFT 默认值为 256 （4） Pxx,w=pmtm( )：输出参数 w 和估计 PSD 的位置一一对应的归一化角频率，单位 rad/sample ，范围如下：若 x 为实信号，则 w 的范围为 0,Pi; 若 x 为复信号，则 w 的范围为 0,2*Pi; 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - -

20、 - - - - - - 第 9 页，共 18 页 - - - - - - - - - 名师精编优秀资料Pxx,f=pmtm( ， Fs): 同时返回和估计PSD 的位置一一对应的线性频率f, 参数 Fs 为采样频率。若 x 为实信号，则 f 的范围为 0,Fs/2; 若 x 为复信号，则 f 的范围为 0,Fs; （5） Pxx,f=pmtm( ， Fs,mehod): Method 有：adapt:Thomson 自适应非线性组合算法，默认值Unity: 相同加权的线性组合Eigen:特征值加权的线性组合（6） Pxx,Pxxc,f=pmtm( ，Fs,mehod,p): P:01，指定

21、PSD 的置信度； Pxxc(:,1)是置信度区间的下部分的数值，Pxxc(:,2)是置信度区间的上部分的数值。（7） Pxx,Pxxc,f=pmtm(x ，dpss_params,NFFT,Fs,mehod,p): （8） =pmtm( ， twoside/oneside )：应用说明：利用Multitaper进行 PSD 估计，并比较 NW 取不同数值时的结果。Fs=1000; nFFT=1024; t=0:1/Fs:1; x=sin(2*pi*200*t)+randn(size(t); P1,f=pmtm(x,2,nFFT,Fs); P2,f=pmtm(x,4,nFFT,Fs); P3,

22、f=pmtm(x,10,nFFT,Fs); Pxx1=10*log10(P1); Pxx2=10*log10(P2); Pxx3=10*log10(P3); subplot(3,1,1); plot(f,Pxx1); xlabel(NW=2); subplot(3,1,2); plot(f,Pxx2); xlabel(NW=4); subplot(3,1,3); plot(f,Pxx3); xlabel(NW=10); 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第

23、 10 页，共 18 页 - - - - - - - - - 名师精编优秀资料050100150200250300350400450500-40-200NW=2050100150200250300350400450500-40-200NW=4050100150200250300350400450500-40-200NW=10可以看出：NW 数值越大，曲线越平滑， 说明方差比较小。 但同时看到波峰变宽，这说明频谱泄露增大了。2、MUSIC 法估计Pmusic 函数：实现 MUSIC 法的功率谱计算调用方式(1) s=pmusic(x,p): 用 MUSIC 法对离散时间信号x 进行功率谱估计。

24、p 是信号 x中包含的复数正弦波信号的个数，如果x 是一个数据矩阵，则对矩阵的每一列都进行功率谱估计。注意：为了返回实信号的全部功率谱值，需要使用另一个输入参数 whole. （3） s=pmusic(r,p, corr ): r 来指定自相关矩阵。 对于实信号而言， 默认情况下Pmusic 返回功率谱估计的一半； 对于复信号， 返回全部的功率谱估计值。（4） s=pmusic(x,p,NFFT): （5） Pxx,w=pmusic( ) （6） Pxx,w=pmusic( ,Fs) 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资

25、料 - - - - - - - - - - - - - - - 第 11 页，共 18 页 - - - - - - - - - 名师精编优秀资料（7） s,f=pmusic( ,NW,noverlap) NW 默认值为 2*p, 参数 noverlap 的 默认值为 NW-1 （8） s,w,v,e=pmusic()： 输出参数 v 为以矩阵，其列是与噪声子空间一一对应的特征值所组成的向量；而e为相关矩阵的特征值向量。应用说明例：用 MUSIC 法进行 PSD估计，结果如图所示：randn(state); n=0:99; s=exp(i*pi/2*n)+2*exp(i*pi/4*n)+exp(

26、i*pi/3*n)+randn(1,100); X=corrmtx(s,7,mod); % 使用改进的协方差方法估计互相关矩阵s,w,e,v=pmusic(X,4,whole); pmusic(X,4,whole); 00.20.40.60.811.21.41.61.8-20-10010203040506070Normalized Frequency ( rad/sample)Power(dB)Pseudospectrum Estimate via MUSIC3、特征向量（ AV）法估计（也是基于矩阵特征分解的一种功率谱估计的非参数方法，它主要适用于混有在噪声的正弦信号的功率谱估计）Peig

27、函数：实现特征向量法的功率谱估计调用方式（1）s=peig(x,p): 用特征向量法对离散时间信号x 进行功率谱估计。 p 是信号x 中包含的复数正弦波信号的个数。(2)s=peig(r,p, corr ) (3)s=peig(x,p,nFFT) (4)Pxx,w=peig( ) (5) Pxx,f=peig( ,Fs) (6)s,f=peig(， NW ，noverlap) 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 12 页，共 18 页 - - - - -

28、 - - - - 名师精编优秀资料(7)s,w,v,e=peig() (8)peig() 应用说明例：用特征向量法进行PSD 估计，结果如下图所示：randn(state,1); n=0:99; s=exp(i*pi/2*n)+2*exp(i*pi/4*n)+exp(i*pi/3*n)+randn(1,100); X=corrmtx(s,7,mod); s,w,e,v=pmusic(X,4,whole); peig(X,4,whole); 00.20.40.60.811.21.41.61.8-30-20-10010203040506070Normalized Frequency ( rad/s

29、ample)Power(dB)Pseudospectrum Estimate via Eigenvector MethodMUSIC算法m=sqrt(-1);delta=0.101043;a1=-0.850848;sample=32; %number of sample spotp=10; %number of sample spot in coef method;名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 13 页，共 18 页 - - - - - - - -

30、 - 名师精编优秀资料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+juinurn= normrnd(0,delta/2,1,sample);uin= normrnd(0,delta/2,1,sample);un=urn+m*uin;%? znfor n=1:sample-1 zn(1)=un(1); zn(n+1)=-a1*zn(n)+un(n+1);end%? xnfor n=1:samplexn(n)=2*co

31、s(2*pi*f1*(n-1)+2*cos(2*pi*f2*(n-1)+2*cos(2*pi*f3*(n-1)+sqrt(2)*real(zn(n);end%*x=xn;for k=0:1:sample-1 s=0;for n=1:sample-k, s=s+conj(x(n,1)*x(n+k,1); %calculate the value of rxx end rxx(1,k+1)=(1/sample)*s;end Rx=zeros(sample,sample); Rx=toeplitz(rxx(1,1:32); U,S,V=svd(Rx); Pmusicf=zeros(1,1/fstep

32、+1); ei=zeros(1,sample);for i=1:length(f)for j=1:sample ei(1,j)=exp(-2*pi*(j-1)*f(i)*m);end ; sum=0;for k=7:sample sum=sum+abs(ei*V(:,k)2;名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 14 页，共 18 页 - - - - - - - - - 名师精编优秀资料end Pmusicf(1,i)=10*log10(1/sum);en

33、d figureplot(f,Pmusicf);-0.5-0.4-0.3-0.2-0.100.10.20.30.40.5-20-15-10-50510名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 15 页，共 18 页 - - - - - - - - - 名师精编优秀资料-90-60-300306090-60-50-40-30-20-100angle (degree)magnitude(dB)derad = pi/180; % deg - rad radeg =

34、180/pi; twpi = 2*pi; kelm = 8; % 阵列数量dd = 0.5; % space d=0:dd:(kelm-1)*dd; % iwave = 3; % number of DOA theta = 10 30 60; % 角度snr = 10; % input SNR (dB) n = 500; % A=exp(-j*twpi*d.*sin(theta*derad);% direction matrix S=randn(iwave,n); X=A*S; X1=awgn(X,snr,measured); Rxx=X1*X1/n; InvS=inv(Rxx); % EV,

35、D=eig(Rxx);% EVA=diag(D); EVA,I=sort(EV A); EVA=fliplr(EV A); EV=fliplr(EV(:,I); % MUSIC for iang = 1:361 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 16 页，共 18 页 - - - - - - - - - 名师精编优秀资料angle(iang)=(iang-181)/2; phim=derad*angle(iang); a=exp(-j*twpi*d*s

36、in(phim).; L=iwave; En=EV(:,L+1:kelm); SP(iang)=(a*a)/(a*En*En*a); end % SP=abs(SP); SPmax=max(SP); SP=10*log10(SP/SPmax); h=plot(angle,SP); set(h,Linewidth,2) xlabel(angle (degree) ylabel(magnitude (dB) axis(-90 90 -60 0) set(gca, XTick,-90:30:90); grid on 2、用 ESPRIT 方法 求 DOA format long N=200;% 快拍

37、数doa=20 40/180*pi; %信号到达角，w=pi/4 pi/3;% 信号频率M=8; %阵元数P=length(w); % 信号个数，也可以用特征分解的大特征值来决定lambda=150;% 波长d=lambda/2;% 阵元间距snr=15;% 信噪比% 导向向量B=zeros(P,M); for k=1:P B(k,:)=exp(-j*2*pi*d*sin(doa(k)/lambda*0:M-1); end B=B; % 导向向量xx=2*exp(j*(w*1:N); x=B*xx; %噪声平均因为 matlab 产生的噪声不太好，不是严格意义上的白噪声，% 平均后结果更好pp

38、,ppp=size(x); 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 17 页，共 18 页 - - - - - - - - - 名师精编优秀资料xx=zeros(pp,ppp); cycle=200; for kkk=1:cycle xx=xx+awgn(x,snr); end x=xx/cycle; %噪声平均结束R=x*x; %数据协方差矩阵%针对相干源的时候进行平衡J=fliplr(eye(M); R=R+J*conj(R)*J; %以下是 ESPRI

39、T 程序Rxx=R(1:M-1,1:M-1);%M-1维的自相关函数Rxy=R(1:M-1,2:M);%M-1维的互相关函数b=zeros(1,M-2);eye(M-2); b=b zeros(M-1,1); Cxx=Rxx-min(eig(Rxx)*eye(M-1); Cxy=Rxy-min(eig(Rxx)*b; a=eig(Cxx,Cxy); %找出最接近 1 的 a 值其对应的角度即为a1=abs(abs(a)-1); for i=1:P c,d=min(a1); a1(d)=1000; bb(i)=a(d); a(d)=1000; end if P1 disp(The angles of signals are) else disp(The angle of signal is) end DOA=asin(angle(bb)/pi)/pi*180 名师归纳总结 精品学习资料 - - - - - - - - - - - - - - -精心整理归纳 精选学习资料 - - - - - - - - - - - - - - - 第 18 页，共 18 页 - - - - - - - - -