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