求助MATLAB中产生高斯色噪声！

peterjiasheng

我想在MATLAB中产生高斯色噪声，首先用randn产生高斯白噪声，然后通过一个低通滤波器，得到高斯色噪声，我想请问下，在MATLAB中如何实现，希望达人能够指教？

飞飞edifier

色噪声？是指有色噪声吗？？

peterjiasheng

就是高斯色噪声呀，

ChaChing

看看这些有没用!?
http://www.mathworks.com/matlabcentral/fileexchange/18269

1. % Function to generate colored noise using Euler-maruyama method
   
2. % For measuring,the accuracy of a numerical simulation to an SDE, the
   
3. % strong order of convergence is used.
   
4. % EM method converges with a strong order =.5
   
5. 
   
6. % Basically to compute the path for color noise we solve the differential
   
7. % equation : dx=axdt+bw(t) where w(t) is wiener process
   
8. % b defines the intensity of white noise
   
9. 
   
10. % Input: a,b (intensity of white noise) and T (the upper limit to time
    
11. % vector), N (the no. of samples which essentially defines the step size )
    
12. % Output: W(t)( Color noise) ,t (time vector)
    
13. 
    
14. % Example colornoise(1,1,1,1000)
    
15. 
    
16. % Author : Abhirup lahiri (abhiruplahiri@yahoo.com)
    
17. 
    
18. function [W,t]= colornoise(a,b,T,N)
    
19. randn('state',100);      %fixing the state of generator
    
20. dt=T/N;
    
21. dW=sqrt(dt)*randn(1,N);
    
22. 
    
23. R=4;                      %fixed for this computation
    
24. L=N/R;
    
25. Xem=zeros(1,L);
    
26. Xzero=0;
    
27. Xtemp=Xzero;
    
28. Dt=R*dt;
    
29. 
    
30. for j=1:L
    
31.     Winc=sum(dW(R*(j-1)+1:R*j));
    
32.     Xtemp=Xtemp+a*Dt*Xtemp+b*Winc;
    
33.     Xem(j)=Xtemp;
    
34. end
    
35. 
    
36. W=[Xzero,Xem];
    
37. t=[0:Dt:T]
    
38. 
    
39. %plot(t,W)
    

复制代码

[ 本帖最后由 ChaChing 于 2009-10-1 08:27 编辑 ]

ChaChing

看看这些有没用!?
http://www.engr.uky.edu/~donohue/ee513/mfiles/MatlabEE513.htm

1. %  This script generates colored noise by filtering white noise with a
   
2. %  Butterworth bandpass filter. Figures are then generated to compare the
   
3. %  PSD of the colored noise to the filter transfer function magnitude.  In
   
4. %  addition, the autocorrelation is taken of the input white noise sequence
   
5. %  and the output colored noise sequence for comparisons.
   
6. %
   
7. %  Written by Kevin D. Donohue (donohue@engr.uky.edu) March 2004
   
8. 
   
9. fs = 8000;  %  Sampling Frequency
   
10. dur = 20;     %  Sound duration in seconds
    
11. ord = 6;
    
12. plen = 32;  % PSD segment length
    
13. %  Bandlimit on the filter
    
14. f1 = [500];  % lower bandlimit
    
15. f2 = [1500]; % corresponding upper bandlimit
    
16. % colors for the plots
    
17. col = ['g', 'r', 'b', 'k', 'c', 'b'];
    
18. no = randn(1,round(fs*dur));  %  Generate noise signal
    
19. [b,a] = butter(ord,2*[f1 f2]/fs);  % Generate filter
    
20. % perform filter operation to get colored noise
    
21. cno = filter(b,a,no);
    
22. % Compute the PSD
    
23. [p, fax] = pwelch(cno,hamming(plen),fix(plen/2),2*plen, fs); % Compute PSD of noise
    
24. figure(1); lh = plot(fax,abs(fs*p/2),col(1)) % Plot PSD
    
25. set(lh,'LineWidth',2)  %  Make line thicker
    
26. hold on
    
27. %  Find filter transfer function
    
28. [h,fq] = freqz(b,a,2*plen,fs);
    
29. plot(fq,abs(h).^2,col(2))
    
30. hold off
    
31. % Label figure
    
32. xlabel('Hertz', 'Fontsize', 14)
    
33. ylabel('PSD', 'Fontsize', 14)
    
34. 
    
35. %  Compute autocorrelation of input sequence
    
36. mxlag = fs*.02;  %  Only compute lags up to 20 milliseconds
    
37. [acwno, lagwno] = xcorr(no, mxlag, 'coef');
    
38. %  Plot lags
    
39. figure(2); plot(1000*lagwno/(fs),acwno)
    
40. xlabel('milliseconds','Fontsize', 14);  ylabel('Correlation coefficient','Fontsize', 14)
    
41. 
    
42. %  Compute autocorrelation of output sequence
    
43. mxlag = fs*.02;  %  Only compute lags up to 20 milliseconds
    
44. [accno, lagcno] = xcorr(cno, mxlag, 'coef');
    
45. %  Plot lags
    
46. figure(3); plot(1000*lagcno/(fs),accno)
    
47. xlabel('milliseconds','Fontsize', 14);  ylabel('Correlation coefficient','Fontsize', 14)
    

复制代码

peterjiasheng

谢谢你！:loveliness: :loveliness:

flexibility

没看明白.........

ChaChing

不见得会懂, 但总需要知道那里没看明白?

guaiguaigui

非常感谢啊

熊星星星

参数a是？？
参数b是指高斯白噪声的信噪比吗？？

熊星星星

Function to generate colored noise using Euler-maruyama method
为什么用不了？Attempt to execute SCRIPT colornoise as a function:

ChaChing

回复 11 # 熊星星星 的帖子

Ref: 常见的程序出错问题整理 http://forum.vibunion.com/thread-46001-1-1.html
   6F->Attempt to execute SCRIPT conv as a function
From http://forum.vibunion.com/home-space-uid-63979-do-blog-id-18250.html