通信人家园

标题: matlab 学习  [查看完整版帖子] [打印本页]

---

时间:  2011-4-22 01:13

作者: timthorpe     标题: matlab 学习

从今天起开始学习MATLAB了


指数序列x=0.8^n*u(n)的离散时间傅立叶变换(1)
w=[0:500]*pi/500;            %  [0,pi]轴等分成501个点
>> X=exp(j*w)./(exp(j*w)-0.2*ones(1,501));
>> magX=abs(X);            angX=angle(X);           %    计算X 的幅度和相位
>> realX=real(X);  imagX=imag(X);                     %   计算X 的实部和虚部
>> subplot(2,2,1); plot(w/pi,magX);  grid
>> xlabel('frequency in pi units'); ylabel('Magnitude'); title('Manitude Part')
>> subplot(2,2,3);plot(w/pi,angX); grid;
>> xlabel('frequncy in pi units'); ylabel('Radians'); title('Angle Part')
>> subplot(2,2,2);plot(w/pi,realX);grid
>> xlabel('frequency in pi units'); ylabel('Real Part'); title('Real Part')
>> subplot(2,2,4);plot(w/pi,imagX);grid
>> xlabel('frequency in pi units');ylabel('Imaginary'); title('Imaginary Part')



[ 本帖最后由 timthorpe 于 2011-4-22 02:18 编辑 ]

---

时间:  2011-4-22 02:18

作者: timthorpe     标题: 指数序列x=0.8^n*u(n)的离散时间傅立叶变换(2)

-------------------------MATLAB脚本-----------------------
        >> n=0:6; x=[4,3,2,1,2,3,4];                        定义频谱计算
                                                              的频率间隔
        >> k=0:500; w=(pi/500)*k;
        >> X=x*(exp(-j*pi/500)).^(n'*k);  //安排为行向量，                        DTFT


        >> magX=abs(X); angX=angle(X);                     分别计算幅度，相
        >> realX=real(X); imagX=imag(X);                    位，实部和虚部
        >> subplot(2,2,1);plot(k/500,magX);grid
        >> xlabel('w');ylabel('幅度');title('Magnitude Part')
        >> subplot(2,2,3);plot(k/500,angX/pi);grid
        >> xlabel('w');ylabel('相位(rad)');title('Angel Part')
        >> subplot(2,2,2);plot(k/500,realX);grid
        >> xlabel('w');ylabel('实部');title('Real Part')
        >> subplot(2,2,4);plot(k/500,imagX);grid
        >> xlabel('w');ylabel('虚部');title('Imaginary Part')

[ 本帖最后由 timthorpe 于 2011-4-22 02:21 编辑 ]

---

时间:  2011-4-22 02:42

作者: timthorpe     标题: 指数序列x=0.8^n*u(n)的离散时间傅立叶变换(3)

设x(n)=0.9exp(jπ/3)n，0≤n≤10,求X(ejw)并研究它的周期性
因为x(n)是复值序列，故只满足周期性质
我们在[ -2π， -2π]之间的401个频率点上求出并画出它的
DTFT以观察它的周期性
>> %x(n)=0.9*exp(j*pi/3)^n,0<n10;
>> n=0:10;x=(0.9*exp(j*pi/3)).^n;
>> k=-200:200;w=(pi/100)*k;
>> X=x*(exp(-j*pi/100)).^(n'*k);
>> magX=abs(X);angX=angle(X);
>> subplot(2,1,1);plot(w/pi,magX);grid
>> subplot(2,1,2);plot(w/pi,angX);grid
>> xlabel('w frequency in units of pi');ylabel('reaisu.pi');title('Angle part');

---

时间:  2011-4-22 02:49

作者: timthorpe     标题: 用MATLAB验证DTFT的性质（以频移性和时移性为例）

用x(n)=cos(πn/2)，0≤n≤100和y(n)=exp(jπn/4) x(n)验证频移性
利用MATLAB分别作出x(n)和y(n)的DTFT的幅度与相位
>> n=0:100;x=cos(pi*n/2);
>> k=-100:100;w=(pi/100)*k;
>> X=x*(exp(-j*pi/100)).^(n'*k);
>> y=exp(j*pi*n/4).*x;
>> Y=y*(exp(-j*pi/100)).^(n'*k);
>> subplot(2,2,1);plot(w/pi,abs(X)); grid; axis([-1,1,0,60]);
>> xlabel('frequency in pi units'); ylabel('|X|'); title('Magnitude of X');
>>subplot(2,2,2);plot(w/pi,angle(X)/pi); grid; axis([-1,1,-1,1]);
>> xlabel('frequency in pi units'); ylabel('radiands/pi'); title('Angle of X');
>> subplot(2,2,3);plot(w/pi,abs(Y)); grid; axis([-1,1,0,60]);
>> xlabel('frequency in pi units'); ylabel('|Y|'); title('Magnitude of Y');
>> subplot(2,2,4);plot(w/pi,angle(Y)/pi); grid; axis([-1,1,-1,1]);
>> xlabel('frequency in pi units'); ylabel('radiands/pi'); title('Angle of Y');

[ 本帖最后由 timthorpe 于 2011-4-22 02:51 编辑 ]

附件: pict.jpg (2011-4-22 02:51, 117.14 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM0MTcwfDJiZGFmYmY5fDE3MzE4Njg0MTl8MHww

---

时间:  2011-4-22 03:06

作者: timthorpe     标题: 验证时移性： 数值法

设x(n)是在0≤n≤10内，[0，1]间均匀分布的随机序列，
并令y(n)=x(n-2)，用MATLAB验证时移性。
>> x=rand(1,11);n=0:10; % 产生随机序列x(n)
>> k=0:500;w=(pi/500)*k; % 定义频率变量w的范围
>> X=x*(exp(-j*pi/500)).^(n'*k); % 计算x(n)的DTFT X(ejw)
>> Y_check=(exp(-j*2).^w).*X; % 用性质计算x(n-2)的DTFT
>> [y,m]=sigshift(x,n,2)； % 产生x(n)的移位序列y(m)
>> Y=y*(exp(-j*pi/500)).^(m'*k); % 用定义求y(m)的DTFT Y(ejw)
>> error=max(abs(Y-Y_check)) % 求两者之间的差的最大值
error =

  8.6483e-015

---

时间:  2011-4-22 19:24

作者: youqiqing

这东西在学校是学过一点  现在都忘记我曾经学过

---

时间:  2011-4-22 19:47

作者: chenaijun     标题: 支持一下

希望分享更多实用的matlab代码，谢谢！

---

时间:  2011-4-22 20:56

作者: shigangl

我是学了等于没学 啊

---

时间:  2011-4-22 21:36

作者: sunnynsj

d

---

时间:  2011-4-22 22:19

作者: timthorpe     标题: 模拟信号傅立叶变换与采样信号的离散时间傅立叶变换

令xa(t)= exp(-1000 |t|),求出并绘制其傅立叶变换Xa(jΩ)。用fs=5kHz进行采样，求出并画出离散时间傅立叶变换X(ejw)。


Dt=0.00005;t=-0.005t:0.005;xa=exp(-1000*abs(t));%模拟信号
Wmax=2*pi*2000;K=500;k=0:1:K;W=k*Wmax/K;
Xa=xa*exp(-j*t'*W)*Dt;Xa=real(Xa);%连续时间傅立叶变换
W=[-fliplr(W),W(2:501)];%频率从-Wmax to Wmax
Xa=[fliplr(Xa),Xa(2:501)];%Xa介于-Wmax和Wmax间
subplot(2,2,1);plot(t*1000,xa);xlabel('时间(毫秒)');
ylabel('xa(t)');title('模拟信号')
subplot(2,2,2);plot(W/(2*pi*1000),Xa*1000);xlabel('频率(kHz)');
ylabel('Xa(jw)');title('连续时间傅立叶变换')
Ts=0.0002;n=-25:1:25;x=exp(-1000*abs(n*Ts));%离散信号
K=500;k=0:1:K;w=pi*k/K;
X=x*exp(-j*n'*w);X=real(X);%离散时间傅立叶变换
w=[-fliplr(w),w(2:K+1)];
X=[fliplr(X),X(2:K+1)];
subplot(2,2,3);stem(n*Ts*1000,x);
xlabel('时间(毫秒)');gtext('Ts=0.2毫秒');
ylabel('x1(n)');title('离散信号')
subplot(2,2,4);plot(w/pi,X);xlabel('频率(弧度)');
ylabel('X1(w)');title('离散傅立叶变换')


[ 本帖最后由 timthorpe 于 2011-4-22 22:28 编辑 ]

附件: sample.jpg (2011-4-22 22:26, 79.21 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM0Mjc0fGQzYTBjNjZifDE3MzE4Njg0MTl8MHww

---

时间:  2011-4-22 22:34

作者: timthorpe     标题: MATLAB 常用函数

[ 本帖最后由 timthorpe 于 2011-4-22 22:35 编辑 ]

---

时间:  2011-5-5 12:59

作者: timthorpe     标题: Tom, Dick, and Mary discover the DFT.pdf

---

时间:  2011-5-5 13:03

作者: timthorpe

Education is the process of telling smaller and smaller lies.


a signal
cannot  be both time-limited,   and
frequency bandlimited,

[ 本帖最后由 timthorpe 于 2011-5-5 13:44 编辑 ]

---

时间:  2011-5-9 14:37

作者: timthorpe     标题: 二阶时域响应曲线

x=0:0.1:20;
y=1-1/sqrt(1-0.3^2)*exp(-0.3*x).*sin(sqrt(1-0.3^2)*x+acos(0.3));
plot(x,y)
h=line(0,0,'color','red','marker','.','markersize',40,'erasemode','xor');
for i=1:length(x)
    set(h,'xdata',x(i),'ydata',y(i));
    pause(0.005)
    drawnow
end

1. %二阶系统时域曲线
   
2. function y=erjiexiangying(zeta)
   
3. %step response of quadratic system
   
4. %zeta 阻尼系数
   
5. %y时域响应
   
6. x=0:0.1:20;
   
7. y=1-1/sqrt(1-zeta^2)*exp(-zeta*x).*sin(sqrt(1-zeta^2)*x+acos(zeta));
   
8. plot(x,y)

复制代码

1. n=20;
   
2. for i=1:n;
   
3.     y=erjiexiangying((i-1)/20);
   
4.     axis([0,20,0,1.5]);
   
5.     M(i)=getframe;
   
6.    
   
7. end;
   
8. movie(M,3)

复制代码

[ 本帖最后由 timthorpe 于 2011-5-9 14:56 编辑 ]

附件: quadratic_reponse.jpg (2011-5-9 14:41, 12.47 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM1NzkzfDM1ODQ4ZWU2fDE3MzE4Njg0MTl8MHww

---

时间:  2011-5-10 15:31

作者: haodiangei     标题: 回复 1# 的帖子

建议，LZ这样效果不是很好。。。加油

---

时间:  2011-5-10 15:45

作者: timthorpe     标题: 连续时间系统零状态响应的求解

1. a = [1 4 4];
   
2. b = [1 3];
   
3. sys = tf(b, a);
   
4. td = 0.01;
   
5. t = 0 : td : 10;
   
6. f = exp(-t);
   
7. y = lsim(sys, f, t);
   
8. plot(t, y);
   
9. xlabel('t(sec)');
   
10. ylabel('y(t)');
    
11. grid on

复制代码

%
sys Transfer function:
    s + 3
-------------
s^2 + 4 s + 4
%计算系统对输入exp(-t)的响应
%

[ 本帖最后由 timthorpe 于 2011-5-10 15:51 编辑 ]

附件: lsim.jpg (2011-5-10 15:47, 32.06 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM2MDM1fGQyNjc4NzAyfDE3MzE4Njg0MTl8MHww

---

时间:  2011-5-10 15:49

作者: timthorpe     标题: 回复 16# 的帖子

能说清楚一点吗 什么效果不好啊?是我学习的方法不好吗 我还是在校生啊，这不过是课上学的，留在这算是复习功课了吧

---

时间:  2011-5-10 15:52

作者: timthorpe     标题: 连续系统冲激响应和阶跃响应求解

---

时间:  2011-5-13 13:51

作者: timthorpe     标题: periodic_prolongation

1. N=24;M=8;m=3;
   
2. n=0:N-1;
   
3. x1=(0.8).^n;
   
4. x2=[(n>=0)&(n<=M)];
   
5. x=x1.*x2;
   
6. xm=zeros(1,N);
   
7. for k=m+1:m+M
   
8.     xm(k)=x(k-m);
   
9. end
   
10. xc=x(mod(n,M)+1);
    
11. xcm=x(mod(n-m,M)+1);
    
12. subplot(4,1,1);stem(n,x,'fill');ylabel('x(n)');
    
13. subplot(4,1,2);stem(n,xm);ylabel('x(n-3)');
    
14. subplot(4,1,3);stem(n,xc);ylabel('x(n)_8');
    
15. subplot(4,1,4);stem(n,xcm);ylabel('x(n-3)_n');

复制代码

[ 本帖最后由 timthorpe 于 2011-5-13 13:52 编辑 ]

附件: periodic_prolongation.jpg (2011-5-13 13:52, 37.37 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM2MzkzfDNmZjU1NjE5fDE3MzE4Njg0MTl8MHww

---

时间:  2011-5-13 14:14

作者: timthorpe     标题: conv_example

1. %x[n]=0.9^n;n=0:19;hn=delta(n-3);
   
2. Nx=20;Nh=10;
   
3. n=0:Nx-1;
   
4. x=(0.9).^n;
   
5. nh=0:Nh-1;
   
6. h=(nh==3);
   
7. ny=0:Nh+Nx-2;
   
8. y=conv(x,h);
   
9. subplot(3,1,1);stem(n,x);
   
10. xlabel('n');ylabel('x(n)');
    
11. subplot(3,1,2);stem(nh,h);
    
12. xlabel('n');ylabel('h(n)');
    
13. subplot(3,1,3);stem(ny,y);
    
14. xlabel('n');ylabel('y(n)');

复制代码

1. %  下面程序是用来实现h和x的卷积得，分别用了filter和conv函数，两者函数得出的结果一样。
   
2. h = [3 2 1 -2 1 0 -4 0 3];       % impulse response
   
3. x = [1 -2 3 -4 3 2 1];            % input sequence
   
4. y = conv(h,x);
   
5. n = 0:14;
   
6. subplot(2,1,1);
   
7. stem(n,y);
   
8. xlabel('Time index n'); ylabel('Amplitude');
   
9. title('Output Obtained by Convolution'); grid;
   
10. x1 = [x zeros(1,8)];
    
11. y1 = filter(h,1,x1);
    
12. subplot(2,1,2);
    
13. stem(n,y1);
    
14. xlabel('Time index n'); ylabel('Amplitude');
    
15. title('Output Generated by Filtering'); grid;

复制代码

[ 本帖最后由 timthorpe 于 2011-5-13 14:22 编辑 ]

附件: conv_example.jpg (2011-5-13 14:15, 27.71 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM2NDAwfDlmNzBkYTU5fDE3MzE4Njg0MTl8MHww

附件: conv_example2.jpg (2011-5-13 14:22, 34.9 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM2NDAyfGI5MDRkNDI4fDE3MzE4Njg0MTl8MHww

---

时间:  2011-5-13 14:33

作者: timthorpe     标题: 连续时间信号卷积

1. function [f,k]=sconv(f1,f2,k1,k2,p)
   
2. %计算连续信号卷积积分f(t)=f1(t)*f2(t)
   
3. %f:    卷积积分f(t)对应的非零样值向量
   
4. %K:    f(t)的对应时间向量
   
5. %f1:   f1(t)的非零样值向量
   
6. %f2:   f2(t)的非零样值向量
   
7. %K1:    序列f1(t)的对应时间向量
   
8. %K2:    序列f2(t)的对应时间向量
   
9. %p:     取样时间间隔
   
10. f1=0.5*(0:0.01:2);f2=0.5*(0:0.01:2);k1=0:0.01:2;k2=0:0.01:2;p=0.01;
    
11. f=conv(f1,f2);                                            %计算序列1与序列2的卷积和
    
12. f=f*p;
    
13. k0=k1(1)+k2(1);                                           %计算序列f非零样值的起点位置
    
14. k3=length(f1)+length(f2)-2;                               %计算卷积和f非零样值得宽度
    
15. k=k0:p:k0+k3*p;                                           %确定卷积和f非零样值的时间向量
    
16. subplot(3,3,1)
    
17. plot(k1,f1)                                               %在子图1绘制f1(t)时域波形图
    
18. title('f1(t)')
    
19. xlabel('t')
    
20. ylabel('f1(t)')
    
21. subplot(3,3,4)
    
22. plot(k2,f2)                                               %在子图2绘制f2(t)时域波形图
    
23. title('f2(t)')
    
24. xlabel('t')
    
25. ylabel('f2(t)')
    
26. subplot(3,3,7)
    
27. plot(k,f);                                                 %画卷积f(t)的时域波形
    
28. h=get(gca,'position');
    
29. h(3)=2.5*h(3);
    
30. set(gca,'position',h)                                     %将第三个子图的横坐标范围扩为原来的2.5倍
    
31. title('      f(t)=f1(t)*f2(t)')
    
32. xlabel('t')
    
33. ylabel('f(t)')

复制代码

[ 本帖最后由 timthorpe 于 2011-5-13 14:34 编辑 ]

附件: sconv.jpg (2011-5-13 14:34, 15.68 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM2NDA3fDk0YTU1MTBmfDE3MzE4Njg0MTl8MHww

---

时间:  2011-5-16 21:02

作者: woshiob

以前学过
久了不用，忘完了。

---

时间:  2011-5-17 14:16

作者: timthorpe     标题: am调制

1. %am调制
   
2. %信源m(t)=sqrt(2)*cos(2*pi*t),载波为s(t)=2*cos(20*pi*t);
   
3. %Signal
   
4. dt=0.01;fmax=1;fc=10;T=5;N=T/dt;
   
5. t=[0:N-1]*dt;
   
6. mt=sqrt(2)*cos(2*pi*fmax*t);%信源
   
7. %AM modulation
   
8. A=2;
   
9. s_am=(A+mt).*cos(2*pi*fc*t);
   
10. %Power Spectrum Density
    
11. [f Xf]=FFT_SHIFT(t,s_am);%调制信号频谱
    
12. PSD=(abs(Xf).^2)/T;%调制信号功率谱密度
    
13. figure(1)
    
14. subplot(211);plot(t,s_am);hold on;%绘制am信号波形
    
15. plot(t,A+mt,'r--');%绘制AM的包络
    
16. title('AM信号及其包络');
    
17. xlabel('t');
    
18. subplot(212);plot(f,PSD);axis([-2*fc 2*fc 0 1.5*max(PSD)]);
    
19. title('AM信号功率谱');xlabel('f');

复制代码

[ 本帖最后由 timthorpe 于 2011-5-17 14:17 编辑 ]

附件: am_modulation.jpg (2011-5-17 14:17, 32.24 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM2ODA0fDk1OTVkOGU4fDE3MzE4Njg0MTl8MHww

---

时间:  2011-5-17 14:34

作者: timthorpe     标题: 抑制载波双边带调制

1. %dsb调制
   
2. %mt=sqrt(2)*cos(2*pi*t),st=cos(2*pi*t)
   
3. %signal
   
4. dt=0.01;
   
5. fmax=1;
   
6. fc=10;
   
7. T=5;
   
8. t=0:dt:T;
   
9. mt=sqrt(2)*cos(2*pi*fmax*t);%信源
   
10. %dsb modulation
    
11. s_dsb=mt.*cos(2*pi*fc*t);
    
12. %power spectrum density
    
13. [f,sf]=FFT_SHIFT(t,s_dsb);
    
14. PSD=(abs(sf).^2)/T;
    
15. %plot DSB AND PSD
    
16. figure (1)
    
17. subplot(211)
    
18. plot(t,s_dsb);hold on
    
19. plot(t,mt,'r--');
    
20. title('DSB 调制信号及其包缝');
    
21. xlabel('t');
    
22. subplot(212);
    
23. plot(f,PSD);
    
24. axis([-2*fc *fc 0 max(PSD)]);
    
25. title('DSB 信号功率谱');
    
26. xlabel('f');

复制代码

附件: dsb_modulation.jpg (2011-5-17 14:35, 85.96 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM2ODA5fDhhZDRmNWNjfDE3MzE4Njg0MTl8MHww

---

时间:  2011-5-18 12:40

作者: timthorpe     标题: am调制与解调

1. %信号的ad调制和解调
   
2. %输入的信号参数
   
3. dt=0.01;
   
4. fm=1;
   
5. fc=10;
   
6. T=5;
   
7. N=T/dt;
   
8. t=[0:N-1]*dt;
   
9. mt=sqrt(2)*cos(2*pi*fm*t); %调制信号 即信号源
   
10. N0=0.01;
    
11. A=2;
    
12. amd=(A+mt).*cos(2*pi*fc*t);%am调制后的信号amt
    
13. B=2*fm;
    
14. noise=GaussNB(fc,B,N0,t);%干扰的高斯噪声
    
15. amd=amd+noise;          % 调制后的信号受到干扰
    
16. figure(1)
    
17. subplot(2,1,1);
    
18. plot(t,amd,t,A+mt,'r--');      %绘出调制后的信号
    
19. 
    
20. legend('调制后信号','包络');
    
21. id_amd=amd.*cos(2*pi*fc*t);%乘以与载波同频同相信号解调
    
22. id_amd=id_amd-mean(id_amd);
    
23. [f,AMF]=FFT_SHIFT(t,id_amd);
    
24. [t,id_amd]=RECT_LPF(f,AMF,B);  %解调后的信号低通滤波
    
25. title('AM 信号');xlabel('t')
    
26. subplot(2,1,2);
    
27. plot(t,id_amd,t,mt/2,'r--');         %绘出解调后的信号 与原信号的对比
    
28. legend('解调信号','原信号');
    
29. title('AM 解调信号');xlabel('t');

复制代码

[ 本帖最后由 timthorpe 于 2011-5-18 13:21 编辑 ]

附件: am_modulation.jpg (2011-5-18 13:21, 42.85 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM2OTM4fGZlMmZiZDMyfDE3MzE4Njg0MTl8MHww

---

时间:  2011-5-20 14:11

作者: timthorpe     标题: 这个很给力啊

1. [x,y,z]=meshgrid(linspace(-1.3,1.3));val=(x.^2 + (9/4)*y.^2 + z.^2 - 1).^3 - x.^2.*z.^3 - (1/9)*y.^2.*z.^3;isosurface(x,y,z,val,0);axis equal;view(-10,24);colormap([1 0.2 0.2])

复制代码

[ 本帖最后由 timthorpe 于 2011-5-20 14:12 编辑 ]

附件: myheart.jpg (2011-5-20 14:12, 14.57 KB) / 下载次数 2
https://www.txrjy.com/forum.php?mod=attachment&aid=MTM3MTcyfDhkNGExMDQyfDE3MzE4Njg0MTl8MHww

---

时间:  2012-3-27 18:29

作者: rushpower

不错

---

通信人家园 (https://www.txrjy.com/)

Powered by C114