matlab与信号处理知识点.docx

matlab与信号处理知识点.docx

《matlab与信号处理知识点.docx》由会员分享，可在线阅读，更多相关《matlab与信号处理知识点.docx（22页珍藏版）》请在冰豆网上搜索。

matlab与信号处理知识点

安装好MATLAB2012后再安装目录下点击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=

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）/N

MATLAB窗函数

矩形窗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.30.61];

a=[10.51];

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=[1234;5678;2478;1023;4567]；

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）行、N^2列的矩阵，该矩阵的列是由矩阵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-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为偶数，则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:

0~1，指定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;%numberofsamplespot

p=10;%numberofsamplespotincoefmethod;

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=normrnd（0,delta/2,1,sample）;

uin=normrnd（0,delta/2,1,sample）;

un=urn+m*uin;

%¼ÆËãzn

forn=1:

sample-1

zn

（1）=un

（1）;

zn（n+1）=-a1*zn（n）+un（n+1）;

end

%¼ÆËãxn

forn=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';

fork=0:

1:

sample-1

s=0;

forn=1:

sample-k,

s=s+conj（x（n,1））*x（n+k,1）;%calculatethevalueofrxx

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）;

fori=1:

length（f）

forj=1:

sample

ei（1,j）=exp（-2*pi*（j-1）*f（i）*m）;

end;

sum=0;

fork=7:

sample

sum=sum+abs（ei*V（:

k））^2;

end

Pmusicf（1,i）=10*log10（1/sum）;

end

figure

plot（f,Pmusicf）;

derad=pi/180;%deg->rad

radeg=180/pi;

twpi=2*pi;

kelm=8;%_阵列数量

dd=0.5;%_space

d=0:

dd:

（kelm-1）*dd;%_

iwave=3;%_numberofDOA

theta=[103060];%_角度

snr=10;%inputSNR（dB）

n=500;%

A=exp（-j*twpi*d.'*sin（theta*derad））;%%%%directionmatrix

S=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））;

%MUSIC_

foriang=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（[-9090-600]）

set（gca,'XTick',[-90:

30:

90]）;

gridon

2、用ESPRIT方法求DOA

formatlong

N=200;%%快拍数

doa=[2040]/180*pi;%%信号到达角，

w=[pi/4pi/3]';%%信号频率

M=8;%%阵元数

P=length（w）;%%信号个数，也可以用特征分解的大特征值来决定

lambda=150;%波长

d=lambda/2;%阵元间距

snr=15;%%信噪比

%%%导向向量

B=zeros（P,M）;

fork=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,ppp]=size（x）;

xx=zeros（pp,ppp）;

cycle=200;

forkkk=1:

cycle

xx=xx+awgn（x,snr）;

end

x=xx/cycle;

%%%%噪声平均结束

R=x*x';%数据协方差矩阵

%针对相干源的时候进行平衡

J=fliplr（eye（M））;

R=R+J*conj（R）*J;

%%以下是ESPRIT程序

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=[bzeros（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）;

fori=1:

P

[c,d]=min（a1）;

a1（d）=1000;

bb（i）=a（d）;

a（d）=1000;

end

ifP>1

disp（'Theanglesofsignalsare'）

else

disp（'Theangleofsignalis'）

end

DOA=asin（angle（bb）/pi）/pi*180