matlab与信号 处理知识点.doc

安装好MATLAB 2012后再安装目录下点击setup.exe会出现 查找安装程序类时出错，查找类时出现异常的错误提示。该错误的解决方法是进入安装目录下的bin文件夹双击matlab.exe对安装程序进行激活。这是可以对该matlab.exe创建桌面快捷方式，以后运行程序是直接双击该快捷方式即可。信号运算1、 信号加MATLAB实现: x=x1+x22、 信号延迟y(n)=x(n-k)3、 信号乘 x=x1.*x24、 信号变化幅度y=k*x5、 信号翻转y=fliplr(x)6、 信号采样和数学描述：y= MATLAB实现： y=sum(x(n1:n2)7、 信号采样积数学描述： MATLAB实现： y=prod(x(n1:n2)8、 信号能量数学描述：MATLAB实现：Ex=sum(abs(x)2)9、 信号功率数学描述：MATLAB实现：Px=sum(abs(x)2)/NMATLAB窗函数矩形窗 w=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发生变化，其旁瓣高度不变。切比雪夫窗：主瓣宽度最小，具有等波纹型，切比雪夫窗在边沿的采样点有尖峰。W=chebwin(n,r)数字滤波器的特性分析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用来指定脉冲信号的采样频率(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个点。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);3、幅频和相频y=abs(x):计算x各元素的绝对值。当x为一个复数时，计算x的复数模。Y=angle(x):计算x向量各元素的复数相位值，单位为弧度。功率谱估（）一、 随机信号处理基础1、 mean函数调用方式：（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)是该矩阵每一个列向量的方差。（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): 估计向量A的自相关函数。C=xcorr(A): 当A为M*N的矩阵时，返回结果为（2M-1）行、N2列的矩阵，该矩阵的列是由矩阵A所有列之间的互相关函数构成。C=xcorr(，scaleopt): 参数scaleopt用来指定相关函数估计所采用的估计方式，即 biased: 有偏估计方式 unbiased: 无偏估计方式 coeff: 对序列进行归一化处理 none： 计算序列的非归一化相关（2） xcov函数：协方差函数估计（3） 相关系数估计计算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为奇数，则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*100*n)+4*sin(2*pi*500*n)+randn(size(n);window=boxcar(length(x);periodogram(x,window,NFFT,Fs);例 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);2、 间接法Fs=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,unbiased);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);Pxx1=10*log10(p1);Pxx2=10*log10(p2);plot(f,Pxx1,r:,f,Pxx2,g-);ylabel(幅值:dB);xlabel(功率谱估计);legend(协方差方法,改进的协方差方法);三、 现在谱估计的非参数方法1、 MTM（Multitaper）法估计时间-带宽的乘积NW，窗口数=2*NW-1Pmtm函数：实现Multitaper法的功率谱估计调用方式（1） Pxx=pmtm(x): 用Multitaper法对离散时间信号x进行功率谱估计，若x为实数，则返回结果为“单边” 功率谱，若x为复数，则返回结果为“双边” 功率谱。（2） Pxx=pmtm(x,NW):（3） Pxx=pmtm(x,NW,NFFT): 参数NFFT用来指定FFT运算所采用的点数。若x为实信号，NFFT为偶数，则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;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，指定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,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); 可以看出：NW数值越大，曲线越平滑，说明方差比较小。但同时看到波峰变宽，这说明频谱泄露增大了。2、MUSIC法估计Pmusic函数：实现MUSIC法的功率谱计算调用方式(1) s=pmusic(x,p): 用MUSIC法对离散时间信号x进行功率谱估计。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)（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(i*pi/3*n)+randn(1,100);X=corrmtx(s,7,mod); %使用改进的协方差方法估计互相关矩阵s,w,e,v=pmusic(X,4,whole);pmusic(X,4,whole);3、特征向量（AV）法估计（也是基于矩阵特征分解的一种功率谱估计的非参数方法，它主要适用于混有在噪声的正弦信号的功率谱估计）Peig函数：实现特征向量法的功率谱估计调用方式（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)(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);MUSIC算法m=sqrt(-1);delta=0.101043;a1=-0.850848;sample=32; %number of sample spotp=10; %number of sample spot in coef method;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 % xn for n=1:sample xn(n)=2*cos(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+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; end Pmusicf(1,i)=10*log10(1/sum); end figureplot(f,Pmusicf);derad = pi/180; % deg - radradeg = 180/pi;twpi = 2*pi;kelm = 8; % 阵列数量dd = 0.5; % space d=0:dd:(kelm-1)*dd; % iwave = 3; % number of DOAtheta = 10 30 60; % 角度snr = 10; % input SNR (dB)n = 500; % A=exp(-j*twpi*d.*sin(theta*derad);% direction matrixS=randn(iwave,n);X=A*S;X1=awgn(X,snr,measured);Rxx=X1*X1/n;InvS=inv(Rxx); %EV,D=eig(Rxx);% EVA=diag(D);EVA,I=sort(EVA);EVA=fliplr(EVA);EV=fliplr(EV(:,I);% MUSICfor iang = 1:361 angle(iang)=(iang-181)/2; phim=derad*angle(iang); a=exp(-j*twpi*d*sin(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 on2、 用ESPRIT方法 求DOAformat longN=200