|
马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
最近 我想做一块关于给定某输入, 以及传递函数的基础上, 计算出结构的响应情况, 我想知道我的程序错误出在哪里? 主要是模拟两个谐频正弦信号x1,x2,再定义两个传递函数H1,H2,然后分别计算a=x1+x2*H1, b=x2+x1*H2 本人很着急,求达人- clc
- clear
- fs=200;
- t=0:1/fs:1;
- N=length(t);
- x1=6*sin(2*pi*10*t)+5*sin(2*pi*20*t)+4*sin(2*pi*30*t)+3*sin(2*pi*40*t)+2*sin(2*pi*50*t)+sin(2*pi*60*t);
- x2=6*sin(2*pi*15*t)+5*sin(2*pi*30*t)+4*sin(2*pi*45*t)+3*sin(2*pi*60*t)+2*sin(2*pi*75*t)+sin(2*pi*90*t);
- % x1=sin(2*pi*10*t)+sin(2*pi*20*t)+sin(2*pi*30*t)+sin(2*pi*40*t)+sin(2*pi*50*t)+sin(2*pi*60*t);
- % x2=sin(2*pi*15*t)+sin(2*pi*30*t)+sin(2*pi*45*t)+sin(2*pi*60*t)+sin(2*pi*75*t)+sin(2*pi*90*t);
- figure
- plot(t,x1)
- hold on
- plot(t,x2,'r')
- X1=fft(x1)/N;X2=fft(x2)/N;
- f=(1:N)*fs/N;
- figure
- plot(f(1:N/2),abs(X1(1:N/2))*2)
- hold on
- plot(f(1:N/2),abs(X2(1:N/2))*2,'r')
- H1=sqrt(1./(300-f.^2).^2+(40.*f).^2).*exp(-j*atan(40.*f./(300-f.^2)));
- H2=sqrt(1./(900-f.^2).^2+(75.*f).^2).*exp(-j*atan(75.*f./(900-f.^2)));
- S1=X2.*conj(X2).*H1./conj(X2);
- for i=1:fix(N/2)
- S1(1,fix(N/2)+i)=0;
- end
- S2=X1.*conj(X1).*H2./conj(X1);
- for i=1:fix(N/2)
- S2(1,fix(N/2)+i)=0;
- end
- AA=2*X1+S1;
- BB=2*X2+S2;
- % AA=X1+X2.*H1;
- % BB=X2+X1.*H2;
- ifftAA=ifft(AA);
- a=real(ifftAA);
- a=a+randn(1,N);
- ifftBB=ifft(BB);
- b=real(ifftBB);
- b=b+randn(1,N);
- figure
- plot(t,a)
- hold on
- plot(t,b,'r')
- AAA=abs(fft(a))*2/N;BBB=abs(fft(b))*2/N;
- figure
- plot(f(1:N/2),AAA(1:N/2))
- hold on
- plot(f(1:N/2),BBB(1:N/2),'r')
复制代码 |
|