求达人，在处理处理传递函数仿真时，出现了问题

sanmuhaifeng

最近  我想做一块关于给定某输入， 以及传递函数的基础上，  计算出结构的响应情况，  我想知道我的程序错误出在哪里？  主要是模拟两个谐频正弦信号x1,x2，再定义两个传递函数H1,H2，然后分别计算a=x1+x2*H1,  b=x2+x1*H2  本人很着急，求达人

1. clc
   
2. clear
   
3. 
   
4. fs=200;
   
5. 
   
6. t=0:1/fs:1;
   
7. N=length(t);
   
8. 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);
   
9. 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);
   
10. 
    
11. % 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);
    
12. % 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);
    
13. 
    
14. figure
    
15. plot(t,x1)
    
16. hold on
    
17. plot(t,x2,'r')
    
18. 
    
19. X1=fft(x1)/N;X2=fft(x2)/N;
    
20. f=(1:N)*fs/N;
    
21. figure
    
22. plot(f(1:N/2),abs(X1(1:N/2))*2)
    
23. hold on
    
24. plot(f(1:N/2),abs(X2(1:N/2))*2,'r')
    
25. 
    
26. H1=sqrt(1./(300-f.^2).^2+(40.*f).^2).*exp(-j*atan(40.*f./(300-f.^2)));
    
27. H2=sqrt(1./(900-f.^2).^2+(75.*f).^2).*exp(-j*atan(75.*f./(900-f.^2)));
    
28. 
    
29. S1=X2.*conj(X2).*H1./conj(X2);
    
30. for i=1:fix(N/2)
    
31.     S1(1,fix(N/2)+i)=0;
    
32. end
    
33. S2=X1.*conj(X1).*H2./conj(X1);
    
34. for i=1:fix(N/2)
    
35.     S2(1,fix(N/2)+i)=0;
    
36. end
    
37. AA=2*X1+S1;
    
38. BB=2*X2+S2;
    
39. 
    
40. 
    
41. % AA=X1+X2.*H1;
    
42. % BB=X2+X1.*H2;
    
43. 
    
44. ifftAA=ifft(AA);
    
45. a=real(ifftAA);  
    
46. a=a+randn(1,N);
    
47. ifftBB=ifft(BB);
    
48. b=real(ifftBB);  
    
49. b=b+randn(1,N);
    
50. 
    
51. 
    
52. figure
    
53. plot(t,a)
    
54. hold on
    
55. plot(t,b,'r')
    
56. 
    
57. AAA=abs(fft(a))*2/N;BBB=abs(fft(b))*2/N;
    
58. figure
    
59. plot(f(1:N/2),AAA(1:N/2))
    
60. hold on
    
61. plot(f(1:N/2),BBB(1:N/2),'r')

复制代码