My goal is to find the spectral density $S_{XX}(f)$ from a time domain noise signal $X(t)$. The approach should work for general time series, but to illustrate what I've been doing I'll keep it simple: we take Gaussian noise with $\mu = 0$ and $\sigma^2 = S_{XX}(f) = 1$, so white noise with a unit spectral density. To show where I get stuck, I'll write down what I've done so far.

I start by generating the Gaussian white noise

UnitData = RandomVariate[NormalDistribution[0, 1], 10^3 + 1];
times = Range[0, 1, 0.001];
Data1 = Table[{times[[i]], UnitData[[i]]}, {i, 1, Length[times]}];

which (using ListPlot) seems to look good.

Okay, well, then we move on to seeing if this indeed has a spectral density of one. My train of thought is as follows.

We have generated a set of $N$ real valued data points $y_r$ (where $r = 1,2,\dots,N$). As we want to analyse the fluctuations of these data points in terms of frequency, we perform a Fourier transformation of our data using Fourier. Mathematically this results in a periodic sequence of points $u_s$ for which $u_{N+s} = u_s$. In addition we know that our input data is real, so that our output is Hermitian: $u_{-s} = u_s^*$. All of the information of our transformation is thus contained in the first $N/2$ Fourier components, as the second half of the components contain the complex conjugates. Finally, in determining the spectral density of our data we are not interested in the phase of the components, only in the absolute values. If we consider the above we can reduce our Fourier components as $N/2$ values of u: \begin{equation} u = \frac{1}{\sqrt{N}}\left(|u_1|,\quad|u_2|+|u_{N}|,\quad|u_3|+|u_{N-1}|,\quad... ,\quad|u_{N/2}|+|u_{N/2+2}|\right) \end{equation} where the factor of $\frac{1}{\sqrt{N}}$ takes the normalization of the data into account, and the addition of the conjugate pairs restores the amplitudes of the signals split into positive and negative frequencies during the Fourier transformation. Note that this splitting does not apply to the first term, which corresponds to $f = 0$.

While we have now transformed our time spectrum data $y_r(t_r)$ to the frequency spectrum data $u_s(f_s)$, we have to find the corresponding frequency array $f_s$. We calculate this from our time data $t_r$ such that \begin{equation} f = \left(0,\,\frac{1}{t_\mathrm{tot}},\,\frac{2}{t_\mathrm{tot}},\,...,\,\frac{1}{2\Delta t}\right) \end{equation} where $\Delta t = t_2 - t_1$ is the measurement time step and $t_{tot} = t_N - t_1$ is the total measurement time.

At this point we have then computed the components $u(f)$ which denote the amplitudes as a function of frequency. In order to obtain the spectral density (or equivalently, the density of the mean square in frequency space $\overline{u_{n}^{2}}(f)$) from the computed $u(f)$ we perform three more steps: first we convert our measurement to the root-mean-square (rms) value by dividing by $\sqrt{2}$, we normalize the spectrum by dividing by the square root of the frequency step size $1/\sqrt{t_{tot}}$ then finally we square the data so that we have \begin{equation} \overline{u_{n}^{2}}(f) = \frac{t_{tot} u(f)^2}{2} \end{equation}

Code wise, I do this with two functions:

DFT[xpts_, ypts_, opts___] := Block[{ans, fpts, freq},
  freq = Range[0, (xpts[[2]] - xpts[[1]])^-1/2, 1/(
    xpts[[-1]] - xpts[[1]])];
  fpts = Fourier[ypts];
  ans = {freq, 1/Sqrt[Length[ypts]] Join[{Abs[fpts[[1]]]},
        Abs[fpts[[i + 1]]] + Abs[fpts[[-i]]], {i, 
         1, (Length[fpts] - 1)/2}]]}\[Transpose];

This guy covers the Fourier transform part, and the next function

NoiseSpectrum[dataname_, dataColor_] := 
  Block[{data, tmax, Δt, stepFreq, maxFreq, datadftResult, 
    SpectralDensityPlot, dataDFT, interpolDFT},
   data = dataname;
   tmax = data[[-1, 1]];
   Δt = data[[2, 1]];
   stepFreq = 1/tmax;
   maxFreq = 1/(2 Δt);
   datadftResult = DFT[Transpose[data][[1]], Transpose[data][[2]]]; 
   interpolDFT = 
    ListInterpolation[Transpose[datadftResult][[2]], {{0, maxFreq}}, 
     InterpolationOrder -> 1];
   dataDFT = Drop[datadftResult, 1]; (*drops DC element*)
   SpectralDensityPlot = ListPlot[
     {#[[1]], (1/(1/Sqrt[tmax]) #[[2]])^2/2} & /@ dataDFT, 
     Frame -> True, 
     FrameLabel -> {"Frequency [Hz]", "Spectral Density, [Arb/Hz]"}, 
     PlotRange -> {{1/tmax, 1/(2 Δt)}, All}, Joined -> True, 
     Joined -> True, PlotStyle -> dataColor]

calculates the neccesary parameters from the input data, and gives a plot of the spectrum. For the above example I use NoiseSpectrum[TestData, Blue] and I get

enter image description here

This is not really what I expected; I thought I'd get a spectrum with a mean of 1. I understand that the fluctuations will always be there of course, but this just seems quite far from what I expected. So my question is, where did I go wrong?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.