基于MATLAB的信号变换.docx
《基于MATLAB的信号变换.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的信号变换.docx（18页珍藏版）》请在冰豆网上搜索。
基于MATLAB的信号变换
将振幅为1的1Hz正弦波和振幅为0.5的5Hz正弦波相加后进行Fourier分析：
>>N=256;dt=0.02;
n=0:
N-1;t=n*dt;
x=sin（2*pi*t）+0.5*sin（2*pi*5*t）;
m=floor（N/1）+1;
a=zeros（1,m）;b=zeros（1,m）;
>>fork=0:
m-1
fori=0:
N-1
a（k+1）=a（k+1）+2/N*x（i+1）*cos（2*pi*k*i/N）;
b（k+1）=b（k+1）+2/N*x（i+1）*sin（2*pi*k*i/N）;
end
c（k+1）=sqrt（a（k+1）.^2+b（k+1）.^2）;
end;
>>subplot（2,1,1）;plot（t,x）;
>>subplot（2,1,2）;plot（（0:
m-1）/（N*dt）,c）;
将上题进行反变换：
>>if（mod（N,2）~=1）
a（m）=a（m）/2;
end;
>>fori=0:
N-1
xx（i+1）=a
（1）/2;
fork=1:
m-1
xx（i+1）=x（i+1）+a（k+1）*cos（2*pi*k*i/N）+b（k+1）*sin（2*pi*k*i/2）;
end;
end;
>>subplot（2,1,1）;plot（（0:
N-1）*dt,x）;
>>subplot（2,1,2）;plot（（0:
N-1）*dt,xx）;
用振幅为0.8的方波进行Fourier分析，并用分析得到的系数求解当k取不同值时的合成信号图：
>>N=200;dt=4/N;
>>forn=1:
N
if（n*dt>=2）
x（n）=0.8;
else
x（n）=-0.8;
end;
end;
>>figure
（1）;
>>subplot（2,1,1）;plot（（1:
N）*dt,x）;holdon;
>>plot（（1:
N）*dt,zeros（1,N）,'r'）;
>>a=zeros（1,N）;
>>b=zeros（1,N）;
>>mm=floor（N/2）+1;
>>fork=0:
mm-1
a（k+1）=0;
b（k+1）=0;
fori=0:
N-1
a（k+1）=a（k+1）+2/N*x（i+1）*cos（2*pi*k*i/N）;
b（k+1）=b（k+1）+2/N*x（i+1）*sin（2*pi*k*i/N）;
end;
c（k+1）=sqrt（a（k+1）.^2+b（k+1）.^2）;
end;
>>subplot（2,1,2）;plot（（0:
mm-1）/（N*dt）,c）;
>>m=input（'输入谐波输入阶数'）;
输入谐波输入阶数5
>>if（m>floor（N/2）+1）
error（'谐波最大阶数必须小于Nyquist频率对应的阶数'）;
end;
>>if（mod（N,2）~=1）
a（mm）=a（mm）/2;
end;
>>fori=0:
N-1
xx（i+1）=a
（1）/2;
fork=1:
m
xx（i+1）=xx（i+1）+a（k+1）*cos（2*pi*k*i/N）+b（k+1）*sin（2*pi*k*i/N）;
end;
end;
>>figure
（2）;
>>plot（（1:
N）*dt,xx,（0:
N-1）*dt,x）;holdon;
>>plot（（1:
N）*dt,zeros（1,N）,'r'）;
Fourier变换程序：
function[Xk]=dfs（xn,N）
n=[0:
1:
N-1];
k=[0:
1:
N-1];
WN=exp（-j*2*pi/N）;
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;
function[xn]=idfs（Xk,N）
n=[0:
1:
N-1];
k=[0:
1:
N-1];
WN=exp（-j*2*pi/N）;
nk=n'*k;
WNnk=WN.^（-nk）;
xn=（Xk*WNnk）/N;
>>xn=[0123];
>>N=4;
>>Xk=dfs（xn,N）
Xk=
6.0000-2.0000+2.0000i-2.0000-0.0000i-2.0000-2.0000i
已知序列x（n）=cos（2*pi*0.24*n）+cos（2*pi*0.26*n）,n=0:
99,试绘制x（n）及其Fourier变换的幅度图：
>>N=100;dt=1;
n=0:
N-1;t=n*dt;
xn=cos（2*pi*0.24*t）+cos（2*pi*0.26*t）;
Xk=dfs（xn,N）;
magXk=abs（Xk）;
phaXk=angle（Xk）;
subplot（2,2,1）;plot（t,xn）;
xx=idfs（Xk,N）;
x=real（xx）;
>>subplot（2,2,2）;plot（t,x）;
>>k=0:
length（magXk）-1;
>>subplot（2,1,2）;plot（k/（N*dt）,magXk*2/N）;
对上式中的X（n）通过补领增长至120个数：
>>N=120;dt=1;
>>n=0:
N-1;t=n*dt;
>>xn=cos（2*pi*0.24*[0:
99]）+cos（2*pi*0.26*[0:
99]）;
>>xn=[xn,zeros（1,N-100）];
>>plot（t,xn）;
Foerier变换的性质：
线性性：
>>N=100;dt=1;
>>n=0:
N-1;t=n*dt;
>>xn1=cos（2*pi*0.24*t）;
>>xn2=cos（2*pi*0.26*t）;
>>Xk1=dfs（xn1,N）;
>>Xk2=dfs（xn2,N）;
>>magXk1=abs（Xk1）;
>>magXk2=abs（Xk2）;
>>phaXk1=angle（Xk1）;
>>phaXk2=angle（Xk2）;
>>k=0:
length（magXk1）-1;
>>subplot（3,1,1）;plot（k/（N*dt）,magXk1*2/N）;
>>k=0:
length（magXk2）-1;
subplot（3,1,2）;plot（k/（N*dt）,magXk2*2/N）;
>>Xk=dfs（xn1+xn2,N）;
>>magXk=abs（Xk）;
>>phaXk=angle（Xk）;
>>k=0:
length（Xk）-1;
>>subplot（3,1,3）;plot（k/（N*dt）,magXk*2/N）;
信号的时移：
functiony=cirshftt（x,m,N）
iflength（x）>N
error（'Nmustbe>=thelengthofx'）;
end
x=[xzeros（1,N-length（x））];
n=[0:
1:
N-1];
n=mod（n-m,N）;
y=x（n+1）;
>>N=128;dt=1;
>>n=0:
N-1;t=n*dt;
>>xn1=cos（0.14*pi*t）;
>>subplot（3,2,1）;plot（t,xn1）;
>>xn2=cirshftt（xn1,20,N）;
>>subplot（3,2,2）;plot（t,xn2）;
>>Xk1=dfs（xn1,N）;
>>magXk1=abs（Xk1）;
>>phaXk1=angle（Xk1）;
>>k=0:
length（magXk1）-1;
>>subplot（3,2,3）;plot（k/（N*dt）,magXk1*2/N）;
>>subplot（3,2,4）;plot（k/（N*dt）,unwrap（phaXk1））;
>>Xk2=dfs（xn2,N）;
>>magXk2=abs（Xk2）;
>>phaXk2=angle（Xk2）;
>>k=0:
length（magXk2）-1;
>>subplot（3,2,5）;plot（k/（N*dt）,magXk2*2/N）;
>>subplot（3,2,6）;plot（k/（N*dt）,unwrap（phaXk2））;
頻移定理：
>>N=100;dt=1;
>>n=0:
N-1;t=n*dt;
>>x=sin（pi*t/2）;
>>k=0:
N-1;f=k/（N*dt）;
>>X=dfs（x,N）;
>>realX=real（X）;
>>imagX=imag（X）;
>>Xrealshft=cirshftt（realX,10,N）;
>>Ximagshft=cirshftt（imagX,10,N）;
>>Y=idfs（Xrealshft+j*Ximagshft,N）;
>>subplot（2,2,1）;plot（t,real（Y））;
>>subplot（2,2,2）;plot（t,imag（Y））;
>>subplot（2,2,3）;plot（t,x.*real（exp（j*2*pi*n*10/N）））;
>>subplot（2,2,4）;plot（t,x.*imag（exp（j*2*pi*n*10/N）））;
奇函数与偶函数与Fourier变换后实部与虚部的关系：
>>N=100;dt=1;n=0:
N-1;t=n*dt;
>>x=cos（2*pi*0.1*t）+0.5*sin（2*pi*0.2*t）;
>>xe=0.5*sin（2*pi*0.2*t）;
>>xo=cos（2*pi*0.1*t）;
>>k=0:
N-1;f=k/（N*dt）;
>>X=dfs（x,N）;
>>XE=dfs（xe,N）;
>>XO=dfs（xo,N）;
>>XR=real（X）;
>>XI=imag（X）;
>>subplot（2,2,1）;plot（t,XR*2/N）;
>>gridon;
>>subplot（2,2,2）;plot（t,XI*2/N）;gridon;
>>subplot（2,2,3）;plot（f,real（XE）*2/N）;gridon;
>>subplot（2,2,4）;plot（f,real（XO）*2/N）;gridon;
褶积定理：
>>N=100;dt=0.01;
>>n=0:
N-1;t=n*dt;
>>f=n/（N*dt）;
>>h=[1zeros（1,N-1）];
>>f0=2;f1=10;
>>x=chirp（t,f0,1,f1）;
>>xh=conv（x,h）;
>>XH=dfs（xh（1:
N）,N）;
>>subplot（2,2,1）;plot（f,real（XH）*2/N）;
>>subplot（2,2,2）;plot（f,imag（XH）*2/N）;
>>X=dfs（x,N）;
>>H=dfs（h,N）;
>>XH1=X.*H;
>>subplot（2,2,3）;plot（f,real（XH1）*2/N）;
>>subplot（2,2,4）;plot（f,imag（XH1）*2/N）;
重叠相加法计算卷积：
functiony=dirconvt（x1,x2,N）
iflength（x1）>N
error（'Nmustbe>=thelengthofx1'）;
end
iflength（x2）>N
error（'Nmustbe>=thelengthofx2'）;
end
x1=[x1zeros（1,N-length（x1））];
x2=[x2zeros（1,N-length（x2））];
m=[0:
1:
N-1];
x2=x2（mod（-m,N）+1）;
H=zeros（N,N）;
forn=1:
1:
N
H（n,:
）=cirshftt（x2,n-1,N）;
end
y=x1*H';
function[y]=ovrlpsav（x,h,N）
Len=length（x）;
M=length（h）;
M1=M-1;L=N-M1;
h=[hzeros（1,N-M）];
x=[zeros（1,M1）,x,zeros（1,N-1）];
K=floor（（Len+M1-1）/（L））;
Y=zeros（K+1,N）;
fori=0:
K
xi=x（i*L+1:
i*L+N）;
Y（i+1,:
）=dirconvt（xi,h,N）;
end
Y=Y（:
M:
N）';
y=（Y（:
））';
>>x1=[1,2,2];x2=[1,2,3,4];
y=dirconvt（x1,x2,5）
x1=[1,2,2];x2=[1,2,3,4];
y=dirconvt（x1,x2,6）
y=
9491414
y=
14914148
采用FFT高速分段卷积的重叠保留法：
function[y]=hsolpsav（x,h,N）
N=2^（ceil（log10（N）/log10
（2）））;
Len=length;
M=length（h）;
M1=M-1;L=N-M1;
h=fft（h,N）;
x=[zeros（1,M1）,zeros（1,N-1）];
K=floor（（Len+M1-1）/（L））;
Y=zeros（K+1,N）;
fori=0:
K
xi=fft（x（i*L+1:
i*L+N））;
Y（i+1,:
）=real（ifft（xi.*h））;
end
Y=Y（:
M:
N）';
y=（Y（:
））';
利用laplace函数：
>>f=sym（'exp（-t）*sin（a*t）'）;
>>L=laplace（f）
L=
a/（（s+1）^2+a^2）
或
>>symsat
L=laplace（exp（-t）*sin（a*t））
L=
a/（（s+1）^2+a^2）
或
>>symsatv
>>f=exp（-t）*sin（a*t）;
>>L=laplace（f,v）
L=
a/（（v+1）^2+a^2）
利用ilaplace函数求F（s）=s^2/（s^2+1）：
>>F=sym（'s^2/（s^2+1）'）;
>>ft=ilaplace（F）
ft=
dirac（t）-sin（t）
或
>>symss
>>ft=ilaplace（s^2/（s^2+1））
ft=
dirac（t）-sin（t）
利用部分展开法求F（s）=（s+2）/（s^3+4*s^2+3*s）：
>>formatrat;
>>B=[12];
>>A=[1430];
>>[r,p]=residue（B,A）
r=
-1/6
-1/2
2/3
p=
-3
-1
0
利用部分展开法求F（s）=（s-2）/（s*（s+1）^2）：
>>B=[1-2];
>>A=conv（conv（[10],[11]）,conv（[11],[11]））;
>>[r,p]=residue（B,A）
r=
2
2
3
-2
p=
-1
-1
-1
0
系统的传递函数Y（s）=（2*s+13）/（s^2+3*s+2），激励信号为x（t）=4×e^（-2）*u（t），起始条件为y（0）=3，y’（0）=4:
>>symsts;
>>Yzis=（2*s+13）/（s^2+3*s+2）;
>>yzi=ilaplace（Yzis）
yzi=
11*exp（-t）-9*exp（-2*t）
>>xt=4*exp（-2*t）*Heaviside（t）;
Warning:
FunctioncallHeavisideinvokesinexactmatchD:
\MATLAB\z\toolbox\symbolic\@sym\heaviside.m.
>>Xs=laplace（xt）;
>>Yzss=Xs/（s^2+3*s+2）;
>>yzs=ilaplace（Yzss）
yzs=
4*（-1-t）*exp（-2*t）+4*exp（-t）
>>yt=simplify（yzi+yzs）
yt=
15*exp（-t）-13*exp（-2*t）-4*t*exp（-2*t）
升级会员 