MATLAB在信号处理中的应用
本帖最后由 zhouyang664 于 2010-12-19 12:02 编辑刚看了网友一个关于信号处理的MATLAB程序,这里奉上我学习过程中总结的一些东东,希望对大家有用!
Ø
基本信号的MATLAB实现:1.
函数:x=;stem(x);%注:若是产生序列可以用stem(x)代替stem(n,x);当然通过改变为值为1时的序列下标,实现 函数的时移;2.
函数:n=-20:20;t=(n>=0);stem(t);%通过改变(n>=i)来实现 函数的时移;3.
单位斜坡函数:n=-20:20;t=n.*(n>0);stem(t);%通过改变n.*(n>i)或(n-i).*(n>i)来实现不同的单位斜坡函数时移;4.
复指数序列:>> n=-10:10;t=0.1+j*0.3;x=exp(t*n);>> subplot(221);stem(real(x));>> subplot(222);stem(imag(x));>> subplot(223);stem(abs(x));>> subplot(224);stem((180/pi*angle(x)));5.
随机序列:使用rand(1,n)和randn(1,n)产生随机序列;6.
MATLAB信号工具箱还提供了一些其他的常用信号,如:SQUARE,SAWTOOTH,SINC,DIRIC,DIRICHLET,RECTPULS和PULSTRAN,具体用法参考help文件。Ø
MATLAB常用函数:real(x):返回复数的实部;imag(x):返回复数的虚部;abs(x):返回复数的模;angle(x):返回复数的相角;rand(1,n):返回长度为n的上均匀分布的随机序列;randn(1,n):返回长度为n的均值为0,方差为1的高斯随机序列,即白噪声序列;Ø
信号的基本运算:1.
信号相加:function =sig_add(x1,n1,x2,n2)% Implements y(n)=x1(n)+x2(n);%x1,x2:序列%n1,n2:序列的起始/终止下标n=min(min(n1),min(n2)):max(max(n1),max(n2));y1=zeros(1,length(n));y2=y1;y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;y=y1+y2;2.
信号相乘:function =sig_mult(x1,n1,x2,n2)% Implements y(n)=x1(n)*x2(n);%x1,x2:序列%n1,n2:序列的起始/终止下标n=min(min(n1),min(n2)):max(max(n1),max(n2));y1=zeros(1,length(n));y2=y1;y1(find((n>=min(n1))&(n<=max(n1))==1))=x1;y2(find((n>=min(n2))&(n<=max(n2))==1))=x2;y=y1.*y2;3.
信号移位:function =sig_shift(x,m,n0)% Implements y(m+k)=x(n)n=m+n0;y=x;4.
序列折叠:function =sig_fold(x,m)% Implements y(n)=x(-n)y=fliplr(x);n=-fliplr(n)5.
序列奇偶性:function=sig_even_odd(x,n)%Real signal decomposition into even and odd partsif(imag(x)~=0)
error('x is not a real sequence');endm=-fliplr(n);m1=min();m2=max();m=m1:m2;nm=n(1)-m(1);n1=1:length(n);x1=zeros(1,length(m));x1(n1+nm)=x;x=x1;xeven=0.5*(x+fliplr(x));xodd=0.5*(x-fliplr(x));下面的用于求奇偶分量的函数的特点在于:它根据输入的序列的中点为对称中心而求奇偶分量:function =circevod(x)if any(imag(x)~=0)
error('The sequence is not real.')
end
N=length(x);
n=0:(N-1);
xeven=0.5*(x+x(sigmod(-n,N)+1));
xodd=0.5*(x-x(sigmod(-n,N)+1));例: n=-5:17;>> x=10*0.8.^n;;>> =circevod(x);>> subplot(311),stem(x);>> subplot(312),stem(xeven);>> subplot(313),stem(xodd);6.
序列圆周移位:function y=cirshift(x,m,N)%x:input sequence%m:shift number%N:show lengthif length(x)>N
error('N must be greater then length(x)');endx=;n=;n=sigmod(n-m,N);y=x(n+1);sigmod函数用来找出周期序列任意位置n所对应的主值有限序列x(n)中的位置m(若周期序列由有限长序列x(n)产生,周期为N)function m=sigmod(n,N)m=rem(n,N);m=m+N;m=rem(m,N);例:由下例可以看出,循环移位不会改变频谱的幅值,只会改变其相位。n=0:10;x=10*0.8.^n;X=DFT(x,11);subplot(321);stem(x);subplot(323);stem(abs(X));subplot(325);stem(angle(X));y=cirshift(x,5,11);Y=DFT(y,11);subplot(322);stem(y);subplot(324);stem(abs(Y));subplot(326);stem(angle(Y));7.
信号的循环/圆周卷积:function y=circonvt(x1,x2,N)%计算循环卷积if length(x1)>N
error('Length(x1) is not greater then N');endif length(x2)>N
error('Length(x2) is not greater then N');endx1=;x2=;m=0:N-1;x2=x2(mod(-m,N)+1);H=zeros(N,N);for n=1:N
H(n,:)=cirshift(x2,n-1,N);endy=x1*H';运行:>> x1=;x2=;>> y=circonvt(x1,x2,5)y =
9
4
9
14
14>> y=circonvt(x1,x2,4)y =
15
12
9
14>> y=circonvt(x1,x2,7)y =
1
4
9
14
14
8
0>> y=circonvt(x1,x2,10)y =
1
4
9
14
14
8
0
0
0
0以上方法是信号序列的时域循环卷积方法,当然也可以使用频域的卷积方法:先将两个序列进行DFT变换,再把DFT变换的结果进行点乘,然后对相乘的结果进行IDFT变换,即得到卷积结果移位: ;MATLAB实现:y=折叠:y(n)=x(-n);MATLAB实现:y=fliplr(x);采样和: ,MATLAB实现:y=sum(x(n1:n2))采样积: ,MATLAB实现:y=prod(x(n1:n2))信号能量: ,MATLAB实现:Ex=sum(abs(x).^2)
因为有些公式是在MathType中编辑的,无法显示,这里附上源文件,参照着学习一下: thanks for sharing!
学习了 谢谢咯。。。 很好的东西{:{39}:} 顺著帖子上浮, 再次学了下!
忘记在那也看过sig_add/sig_mult这两个函数了, 细看下个人感觉好像可改进
在找序列的起始/终止下标, 使用min(n1), max(n1)...等, 好像根本多餘! min(n1)不就是n1(1), max(n1)不就是n1(end)! 难道时间序列可以不是单调递增吗?
所以练习下, 依个人习惯并无特别优化function =sigadd(f1,n1,f2,n2)
n=min(n1(1),n2(1)):max(n1(end),n2(end)); f=zeros(1,length(n));
ppp=n1(1)<=n & n<=n1(end); f(ppp)=f1;
ppp=n2(1)<=n & n<=n2(end); f(ppp)=f(ppp)+f2; {:{51}:} 学习了。
页:
[1]