2
$\begingroup$

I need to generate random numbers with a given PSD. To do so, I found this recipe. In case it's behind a paywall, here's now it works:

Given your power spectrum $P_k$ with $0\leq k<N$ points generated using a timeseries with timestep $h$, you define

$$A_k=\sqrt{\frac{P_k}{2Nh}}\exp{(ia_k)}$$

where $a_k$ are uniformly randomly chosen in $[0,2\pi)$. The desired timeseries is then the Fourier transform of $A_k$. The idea here is that since the PSD throws away the phase information, we simply have to generate random phases for the transform to get a possible realization of the timeseries that generated our PSD.

However, to ensure that all of the resulting time points are real, the author gives the condition that $A_{N-k}=A^*_k$, or equivalently $a_{N-k}=-a_k$.

When I implement this, however, I don't get a real signal, so I seem to be doing something wrong. I see one problem right away: for $k=0$, the FFT assumes a periodic signal, so the only way to satisfy the condition is if $a_0=0$. Then you can generate the rest of the phases as follows (in python):

def gen_phase(N):
    phase = np.zeros(N)
    randoms = np.random.random(size=N/2)*2*np.pi
    phase[1:N/2+1] = randoms
    phase[N/2+1:] = np.flipud(-randoms)
    return phase

The other problem is that it is not possible to satisfy the given condition is $N$ is even, as far as I can tell. This is a pain practically, because the fft is very slow for odd numbers (particularly ones of the form $2^n+1$...)

Can someone suggest what I might be doing wrong? The signal I get looks more or less right ion terms of the amplitude, but it has roughly equivalent real and imaginary parts, so something is wrong. The rest of the code is too simple to be the culprit, I think:

def main():
    psd = read_psd()
    N = len(psd)
    phase = gen_phase(N)
    h = 0.24e-6
    A = np.sqrt(psd/(2*N*h)) * np.exp(1j*phase)
    timeseries = np.fft.fft(A)
    pl.plot(np.real(timeseries))
    pl.show()

if __name__=='__main__':
    main()

EDIT: I misinterepreted. Because $A_k$ depends on $P_k$, I don't have the freedom to choose $A_{N-k}=A^*_k$ in the first place, because I don't have control over the amplitude. Any ideas?

$\endgroup$
  • $\begingroup$ How to translate this averages to time scale? $\endgroup$ – Sabina Feb 18 at 15:54
  • $\begingroup$ Can someone explain the relation between "averages" and time scale for generated time-domain signal from above code? $\endgroup$ – saab Feb 18 at 18:31
1
$\begingroup$

Hopefully this octave implementation will help you:

N = 1e3;
P = ones(N,1);
Ts = 1/1000;
A = sqrt(P/(2*N*Ts)) .* exp( 1j * rand(N,1) );

% DC sample must be real.
n0 = floor( N/2 ) + 1;
A(n0) = P(n0);

% For even-length cases, the first sample has no
% equivalent positive frequency, so we force it
% to be real too.
A(1) = abs( A(1) );

% Force the spectrum to be conjugate symetric.
if( mod(N,2) == 0 )
  A(2:n0-1) = flipud( conj(A(n0+1:end)) );
else
  A(1:n0-1) = flipud( conj(A(n0+1:end)) );
end

% Time domain data.
x = ifft( ifftshift( A ) );
$\endgroup$
  • $\begingroup$ Correct me if I'm wrong, but it looks like you are just dropping half of your power spectrum in this case when you do the conjugate symmetry. Is that accurate? $\endgroup$ – KBriggs Mar 7 '18 at 17:15
  • $\begingroup$ @KBriggs I am, but for the power spectrum to produce a real-valued signal in the time domain, the magnitudes of the dropped values of A must be the equal to the ones that replace them. If P is not symmetric to begin with, then the process will not work. Remember that for a real-valued signal, half of the spectrum is redundant. $\endgroup$ – AnonSubmitter85 Mar 7 '18 at 17:21
  • $\begingroup$ @KBriggs Note that I added something in an edit. You need to force the first sample of A to be real too. Otherwise, its phase carries over to the time-domain for even N. $\endgroup$ – AnonSubmitter85 Mar 7 '18 at 17:29
1
$\begingroup$

You are reversing your phases okay. You should be able to use an even value of N if you set the phase for N/2 to zero.

The problem you have is that your PSD spans the entire N. Unless $ P_k = P_{N-k} $ you are not generating a DFT with complex conjugate symmetry which is what you need to produce a real valued signal.

Also, you should be using the "ifft" instead of the "fft" to do an inverse FFT.

Hope this helps,

Ced

Edit: I didn't see your edit before I posted. It seems you figured this out.

======================

Followup:

For even N:

def gen_phase(N):
    phase = np.zeros(N)
    randoms = np.random.random(size=N/2-1)*2*np.pi
    phase[1:N/2] = randoms
    phase[N/2+1:] = np.flipud(-randoms)
    return phase
$\endgroup$
  • $\begingroup$ But I don't have control over that, do I? My PSD is just a given. Or can I simply duplicate it and use a full PSD length 2N? $\endgroup$ – KBriggs Mar 7 '18 at 17:09
  • $\begingroup$ @KBriggs, No, that won't work because you will get a sequence that is twice as long. The PSD would have to be mirrored, like you did with the phases, and thus be a different PSD. There is a slick trick luring here, I can sense it. Let me give this some more thought. $\endgroup$ – Cedron Dawg Mar 7 '18 at 17:16
  • $\begingroup$ I think I have an idea: my psd actually already throws out half of the frequency information because it is a real signal and so the negative frequency component is redundant. If I just mirror the PSD about 0, I think I can recover it. I'm running it now, but it will take a while because the FFT doesn't do well with odd numbers. $\endgroup$ – KBriggs Mar 7 '18 at 17:51
  • $\begingroup$ @KBriggs, I think you just put your finger on it. For a real valued signal, the PSD has to be symmetric about 0. So if your given PSD isn't, then it can't be the PSD of a real signal. You should be able to test your code with a much smaller N. See my followup for even N. $\endgroup$ – Cedron Dawg Mar 7 '18 at 18:06
  • $\begingroup$ Progress, but still not quite working. The power spectrum comes out more or less correct, but it is still split between real and imaginary parts. This probably comes down to an off-by-one index somewhere. Fun stuff. $\endgroup$ – KBriggs Mar 7 '18 at 19:40
1
$\begingroup$

Thanks to discussion with the other answerers I was able to figure it out. I am posting my own answer here since none of the others have a complete picture, though they certainly helped.

The trick I was missing is that because my input PSD is only defined for positive frequencies, I had to first duplicate the PSD symmetrically about 0. Then, because of the particular ordering of frequency components internal to numpy's fft implementation, an ifftshift was needed to bring it all into line.

Here is the complete code, which works assuming you have your PSD in a .csv file defined for frequencies in [0, fmax].

import numpy as np
from scipy.fftpack import fft, fftshift, ifftshift
from scipy.signal import welch
import pandas as pd
import matplotlib.pyplot as pl

def read_psd():
    psd = pd.read_csv('psd-example.csv',header=None, names=['f','P','I'])
    return psd['f'],psd['P']

def gen_phase(N):
    if N%2==0:
        raise ValueError('THis methofd requires a symmetric PSD, which can only be achieved with an odd number of data points')
    else:
        phase = np.zeros(N)
        randoms = np.random.random(size=N/2)*2*np.pi-np.pi
        phase[:N/2] = randoms
        phase[N/2+1:] = np.flipud(-randoms)
    return phase


def symmetrize_psd(psd):
    N = len(psd)
    symmetric = np.zeros(2*N-1)
    symmetric[0:N] = np.flipud(psd)
    symmetric[N:] = psd[1:]
    return symmetric

def scale_psd(psd, dt):
    N = len(psd)
    phase = gen_phase(N)
    A = np.sqrt(psd/(2*N*dt)) * np.exp(1j*phase)
    return ifftshift(A)

def get_timeseries(A):
    timeseries = np.fft.fft(A)
    return np.real(timeseries)

def main():
    averages = 200
    fi, psd = read_psd()
    psd = psd[:-1]
    fi=fi[:-1]

    sympsd = symmetrize_psd(psd)
    N = len(sympsd)
    timeseries = np.empty(N*averages)
    for i in range(averages):
        A = scale_psd(sympsd, 0.24e-6)
        timeseries[i*N:(i+1)*N] = get_timeseries(A)
    f, pxxr = welch(np.real(timeseries), fs=4166666.67, nperseg=2**np.floor(np.log2(N)+1))
    pl.loglog(f, pxxr)
    pl.loglog(fi, psd)
    pl.ylim((1e-5,100))
    pl.show()

if __name__=='__main__':
    main()
$\endgroup$
  • 1
    $\begingroup$ I don't see why this requires an odd length FFT. $\endgroup$ – AnonSubmitter85 Mar 8 '18 at 4:50
  • $\begingroup$ If I have a PSD defined on [0,fmax] with n points, then mirroring it about 0 requires another n-1 points (every point except 0) for a total of 2n+1 points in the symmetric PSD. Though now that I think about it it could be an even number of points if the 0-frequency bin was not defined, which is a case not covered by my code above. I will add that when time permits. Thanks for the correction! $\endgroup$ – KBriggs Mar 8 '18 at 14:52

Not the answer you're looking for? Browse other questions tagged or ask your own question.