Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
generating noise with given psd

Subject: generating noise with given psd

From: Robert

Date: 8 Dec, 2011 16:09:09

Message: 1 of 4

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?


-------------------

function [time,x]=makenoise(P,dxi)

% Default bin size is 1
if nargin==1
    dxi=1;
end

time=0;
x=0;
XX=0;

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

phase=randn(1,numel(xis))*2*pi; % random phases
X=sqrt(P).*exp(sqrt(-1)*phase); % Generate random spectrum
XX=[X conj(X(end:-1:1))]; % *** This line is necessary to make noise real
x=ifft(XX,Nxis,'symmetric'); % Convert noise to time domain


% Remove bias, normalize variance
x=x-mean(x);
x=x*sqrt(varx/var(x));


% Build time vector
dt=1/ximax;
T=1/dxi;
time=0:dt:T;

end

Subject: generating noise with given psd

From: Greg Heath

Date: 9 Dec, 2011 00:35:55

Message: 2 of 4

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

Subject: generating noise with given psd

From: Robert

Date: 21 Dec, 2011 18:44:08

Message: 3 of 4

Thanks. As far as I can tell, my mistakes were calling randn versus rand to generate random phases (this was just a typo), and failing to scale the integral of the psd by 2 to get the variance (since it's one-sided), and not being careful about whether there were an odd or even number of frequencies.

Subject: generating noise with given psd

From: John

Date: 4 Jan, 2012 19:41:08

Message: 4 of 4

I tried to follow the corrections, but appear to be doing something incorrectly.
The following code will plot a PSD, interpolate the PSD to even spacing, calculate the corresponding time history, then calculate the PSD of the time history and plots the time history and the original and constructred PSD. The constructed PSD does not match the original PSD. Can someone indicate where the error is?

function [time,x]=makenoise(PSD,df)

if nargin == 0
    f_hz_s_m2_hz_20_6_m_s = [0 6.289308176101E0
    5.403450004557E-3 8.647798742138E0
    1.010504967642E-2 1.100628930818E1
    1.626959711968E-2 9.591194968553E0
    2.062368972746E-2 7.861635220126E0
    2.461603317838E-2 5.974842767296E0
    3.150407893538E-2 4.874213836478E0
    3.403974113572E-2 5.031446540881E0
    3.766919606235E-2 3.301886792453E0
    4.382291951509E-2 4.874213836478E0
    4.382064078024E-2 5.503144654088E0
    6.592607784158E-2 1.044025157233E2
    7.137225412451E-2 1.012578616352E2
    8.893674232066E-2 5.345911949686E1
    9.949924801750E-2 3.820754716981E1
    1.049368790448E-1 3.742138364780E1
    1.118397365783E-1 3.223270440252E1
    1.158525886428E-1 2.468553459119E1
    1.216650715523E-1 2.044025157233E1
    1.281862409990E-1 2.059748427673E1
    1.557572235895E-1 1.100628930818E1
    1.662650396500E-1 1.084905660377E1
    1.764156640233E-1 9.276729559748E0
    1.876498268162E-1 8.647798742138E0
    1.992520052867E-1 6.446540880503E0
    2.112096663932E-1 6.132075471698E0
    2.271613800018E-1 3.459119496855E0
    2.481775818066E-1 2.987421383648E0
    2.565086364051E-1 3.616352201258E0
    2.662940935193E-1 2.830188679245E0
    2.728186810683E-1 2.044025157233E0
    2.807891258773E-1 2.201257861635E0
    2.836899553368E-1 1.572327044025E0
    2.894887658372E-1 1.100628930818E0
    2.945583811868E-1 1.886792452830E0
    2.996285662200E-1 2.515723270440E0];
psd_data = sortrows(f_hz_s_m2_hz_20_6_m_s);

% Linearly interpolate at even spacing
% df_hz = mean(diff(psd_data(:,1)));
df_hz = 0.01;
f_hz_min = ceil(min(psd_data(:,1))/df_hz)*df_hz;
f_hz_max = floor(max(psd_data(:,1))/df_hz)*df_hz;
f_hz = [f_hz_min:df_hz:f_hz_max]';
psd_interp = interp1(psd_data(:,1),psd_data(:,2),f_hz);

df_rad_sec = df_hz*2*pi;
[time,y]=makenoise(psd_interp,df_rad_sec);

close all;
figure;
subplot(2,1,1);
plot(time,y,'k-');
xlabel('T'); ylabel('Wave height (m)');
Y = fft(y);
N = length(Y);
absY = abs(Y);
Pyy = absY.^2/N;
Fs=df_rad_sec*N;
PSD = Pyy/Fs;
subplot(2,1,2);
plot(psd_data(:,1),psd_data(:,2),'bo');
hold on;
plot([0:(N-1)]*df_rad_sec/(2*pi),PSD,'k-');
xlabel('Frequency (Hz)');
ylabel('Wave Spectral Density (m^2/Hz)');
legend('20.6 m/s', 'Reconstructred', 'location', 'northeast');
set(gca,'xlim',[0 .3],'XMinorTick', 'on','XTick',[0:.05:.3], ...
    'ylim',[0 110],'YMinorTick', 'on','YTick',[0:10:110], ...
    'yticklabel',{'0','','20','','40','','60','','80','','100',''});

return;


end


% Default bin size is 1
N = length(PSD);
if nargin==1
    % df=1;
    df = 1/N;
end
Fs = df*N;

ximax=(N-1)*df; % Maximum frequency
f=0:df:(ximax); % vector of frequencies
% Nxis=numel(xis); % number of frequencies

varx=sum(PSD*df); % variance of x, calculated from P

phase=rand(size(PSD))*2*pi; % random phases
% X=sqrt(P).*exp(sqrt(-1)*phase); % Generate random spectrum
X = sqrt(N*Fs*PSD).*exp(i*phase);
isEven = abs(mod(N,2)-1);
if isEven
XX=[X; conj(X(end:-1:1))]; % *** This line is necessary to make noise real
else
XX=[X;0; conj(X(end:-1:1))]; % *** This line is necessary to make noise real
end
x=ifft(XX,N,'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
h=1/ximax;
T=1/df;
time=[0:h:T]';

end
 

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us