# Applying Fourier Shift Theorem using Power Spectrum

I have a power spectrum $$P(\omega) = 10/\omega^2$$ and I want to simulate a random time-series realization of this process $$f(t)$$ and a time-shifted copy of that realization $$f(t- \tau)$$. To do this, I want to the Fourier Shift theorem $$\mathcal{F}[f(t - \tau)](\omega) = e^{-2\pi j \omega \tau}\mathcal{F}[f(t)](\omega)$$

To get $$f(t)$$ I first discretize the frequency domain, get the amplitude of the frequency components by taking $$A(\omega) = \sqrt{2P(\omega)\Delta \omega}$$ and then multiply them by a random complex phase uniformly spaced on $$[0, 2\pi]$$. Finally I perform the Fourier Transform and take the real part to get a random realization of $$f(t)$$. This overall looks something like $$f(m\Delta t) = Re\left\{\sum_{n=1}^N A(\omega_n)e^{j n m 2\pi/N} e^{j\theta}\right\}$$ where $$\theta = \mathcal{U}[0,2\pi]$$.

To get the shifted copy, I use the shift theorem and get

$$f(m\Delta t - \tau) = Re\left\{\sum_{n=1}^N A(\omega_n)e^{j n m 2\pi/N} e^{j\theta}e^{-2\pi \omega_n \tau}\right\}$$ where $$\theta = \mathcal{U}[0,2\pi]$$.

P[w_] := 10^1 w^-2
tot = 2^10;
wmax = .1;
dw = wmax/tot;
dt = N[(2 \[Pi])/(2 dw tot)]
An = Table[2 Sqrt[P[n dw] dw], {n, 1, tot}];
theta = RandomVariate[UniformDistribution[{0, 2 \[Pi]}], tot];
Bn = Exp[I theta];
Cn = Table[Exp[-2 \[Pi] I n dw dt 20], {n, 1, tot}];

f1 = Re[Fourier[An*Bn, FourierParameters -> {1, -1}]];
f2 = Re[Fourier[An*Bn*Cn, FourierParameters -> {1, -1}]];

ListPlot[{f1, RotateRight[f2, 20]}, Joined -> True]



Below is the output of this calculation, and not only do the time-series not line up after I've shifted them, but there are also differences in the amplitude. I would like help in understanding why there are differences in amplitude, why the shifts don't line up, and how to correctly apply the Fourier Shift Theorem.

This is all surprisingly tricky, so I show what I would do, rather than attempting to fix your version.

A number of general points:

• You have to use consistent times and frequencies. Here I have unit spacing in time and consequently a total frequency range of 1. Note that the frequencies appear in the spectrum in a shifted order.
• You (may) need to worry about aliasing. Until you know what you are doing, ensure that your spectra drop to zero rapidly at the frequency discontinuity (where the minimum and maximum frequency meet)

Define the parameters of the problem, and a power spectrum that drops off rapidly.

npoints = 256;
bandwidth = 1/24;
spectrum = Exp[-π (f/bandwidth/2)^2];


Define a function that moves the zero frequency to the edge of the data set (or back again, assuming npoints is even).

fftshift = RotateLeft[#, Length[#]/2] &;


It is important that we define the frequency of our spectral values correctly

framp = fftshift[Range[-npoints/2, npoints/2 - 1]/npoints];


Use complex Gaussian noise

noise = {1, I} . RandomReal[NormalDistribution[0, 1], {2, npoints}];


and shape its spectrum

spec1 = noise*Sqrt[spectrum /. f -> framp];

ListLinePlot[Transpose[{fftshift[framp], fftshift[Abs[spec1]^2]}],
PlotRange -> All]


Define a signal with this power spectrum

signal1 = InverseFourier[spec1, FourierParameters -> {1, -1}];

ListLinePlot[Transpose[ReIm[signal1]]]


(You could check its autocorrelation function, but I've omitted that step).

Define a delay

dt = 7.5;


and apply this delay in the frequency domain

phaseterm = Exp[I*2*Pi*framp*dt];
spec2 = spec1*phaseterm;


Create a delayed version of the signal

signal2 = InverseFourier[spec2, FourierParameters -> {1, -1}];


and see that (the real part) is indeed delayed

ListLinePlot[Re /@ {signal1, signal2}]


Zoom in to verify that the delay has the requested value

ListLinePlot[Re /@ {signal1, signal2}, PlotRange -> {{150, 180}, All}]