Path: news.mathworks.com!newsfeed-00.mathworks.com!newsfeed2.dallas1.level3.net!news.level3.com!postnews.google.com!d10g2000vbf.googlegroups.com!not-for-mail From: Greg Heath <heath@alumni.brown.edu> Newsgroups: comp.soft-sys.matlab Subject: Re: generating noise with given psd Date: Thu, 8 Dec 2011 16:35:55 -0800 (PST) Organization: http://groups.google.com Lines: 105 Message-ID: <c22aa5b9-e9af-4fab-a4a4-29ea62a14101@d10g2000vbf.googlegroups.com> References: <jbqnf5$qo$1@newscl01ah.mathworks.com> NNTP-Posting-Host: 68.83.182.65 Mime-Version: 1.0 Content-Type: text/plain; charset=ISO-8859-1 X-Trace: posting.google.com 1323390956 12984 127.0.0.1 (9 Dec 2011 00:35:56 GMT) X-Complaints-To: groups-abuse@google.com NNTP-Posting-Date: Fri, 9 Dec 2011 00:35:56 +0000 (UTC) Complaints-To: groups-abuse@google.com Injection-Info: d10g2000vbf.googlegroups.com; posting-host=68.83.182.65; posting-account=mUealwkAAACvQrLWvunjg50tRAnsNtJR User-Agent: G2/1.0 X-Google-Web-Client: true X-Google-Header-Order: ARLEUHNKC X-HTTP-UserAgent: Mozilla/4.0 (compatible; MSIE 7.0; Windows NT 5.1; AT&T CSM6.0; AT&T CSM 6; GTB7.2; (R1 1.3); (R1 1.6); .NET CLR 1.0.3705; .NET CLR 2.0.50727; .NET CLR 3.0.4506.2152; .NET CLR 3.5.30729),gzip(gfe) Xref: news.mathworks.com comp.soft-sys.matlab:751897 On Dec 8, 11:09 am, "Robert " <robe...@pitt.edu> wrote: > Hi. I have a power spectrum and I'd like to generate (real-valued) noise that has that spectrum, but I'm having trouble getting what I want. > > I attached what I wrote so far, but it adds a high frequency component to the noise, due to the line marked with a ***, which is necessary to make the noise real-valued. Any suggestions? > ------------------- Is the Power Spectrum one or two-sided? What is the Phase Spectrum? What are the probability distributions for Amplitude and Phase ? My cheat-sheet: N = length(x) absX = abs(fft(x)); Pxx = absX.^2/N; % Mean-Squared Power Spectrum PSD = Pxx/Fs % Power Spectral Density % Parceval (df = Fs/N): % Pav = average power Pav = mean(x.^2) = df*sum(PSD) Pav = mean(x.^2) = mean(Pxx) Pav = mean(x.^2) = var(x)+mean(x)^2 > function [time,x]=makenoise(P,dxi) dxi = df % (Less confusing) Is P = Pxx or PSD? Is P single-sided or double-sided? > % Default bin size is 1 > if nargin==1 > dxi=1; > end Why not Fs = 1 ==> df = 1/N ? > time=0; > x=0; > XX=0; Unnecessary. No loops involved. > ximax=(numel(P)-1)*dxi; % Maximum frequency > xis=0:dxi:(ximax); % vector of frequencies > Nxis=numel(xis); % number of frequencies > > varx=sum(P*dxi); % variance of x, calculated from P Only if 1. P = PSD 2. P is double-sided 3. mean(x) = 0 ( See my cheat-sheet). > phase=randn(1,numel(xis))*2*pi; % random phases Phase is odd and periodic. How can it be Gaussian? 1. Use rand 2. Phase should be asymmetric about Nyquist 3. Phase of DC and Nyquist are zero. > X=sqrt(P).*exp(sqrt(-1)*phase); % Generate random spectrum X is not a spectrum X = sqrt(N*Fs*PSD).*exp(i*phase); % ( See my cheat-sheet). > XX=[X conj(X(end:-1:1))]; % *** This line is necessary to make noise real Only if N is odd . Modify for even N. If PSD is one-sided, some of earlier formulas need to be changed... I'll leave that for you. > x=ifft(XX,Nxis,'symmetric'); % Convert noise to time domain check = max(abs(imag(x))) x = real(x); > % Remove bias, normalize variance > x=x-mean(x); > x=x*sqrt(varx/var(x)); > > % Build time vector > dt=1/ximax; Not if PSD is one-sided > T=1/dxi; Not if PSD is one-sided > time=0:dt:T; max(t) = T-dt > end Hope this helps. Greg