0

I'm trying to learn more about noise, power spectral denisty (PSD) and statistical variances. With regards this I'm trying to compute the Power Spectral density of white noise, however, when I do I get a very odd symmetry. my spectrum seems to be symmetric around the central frequency value, which is obviously incorrect. I'm new to using autocorrelations and Power spectral densities so I'd appreciate if someone could help nudge me in the direction of the error.

Code to calculate PSD:

import numpy as np 
from math import sin, pi, log10
from allan_noise import white,pink,brown, violet
import acor
import numpy as np


#METHOD ONE: AUTOCORRELATION FUNCTION METHOD
def psd1(time_dats):
    #auto_corr = acor.function(time_dats)
    auto_corr = np.correlate(time_dats,time_dats,mode="same")

    #To check autocorrelation
    for i in range(len(auto_corr)):
        print i*.0000001 ,auto_corr[i]

    return fft(auto_corr) 

#DEFINE VARIABLES
t = .0001       #simulation time
N = 100000  #number of data points 
dt = t/N    #time interval between data points

#CREATE SIGNAL
sig = white(N)
df = 1.0/(t)
freq_axis = np.arange(0,N/t,df)

spec = psd1(sig)

#OPEN UP OUTPUT FILE
f = open('data/psdtest_f','w')
g = open('data/psdtest_t','w')

#PRINT OUT DATA
for i in range(N):
    f.write('%g %g\n' %(freq_axis[i],log10(abs(spec[i]))))
    g.write('%g %g\n' %(i*dt, sig[i]))

Using this code I produce the following graphs, which can be accessed here https://drive.google.com/#folders/0B8BdiIhqTYw7Vk1EaDFMQW84RHM :

  1. Temporal Profile of the noise before calculations

  2. Autocorrelation function calculated from the temporal profile (I'm aware the x axis scale is wrong but that doesn't contribute to the code any where else

  3. Power Spectral Denisty, beautifully symmetric though not supposed to be

Any help anyone could provide as to whats causing this symmetry would be most helpful. i've tested the code with a simple sin wave signal and I get expected results (with no symmetry)

Leavenotrace
  • 3
  • 1
  • 3
  • This question appears to be off-topic because it is about interpretation of results, not code implementation – jonrsharpe Jun 19 '14 at 11:15
  • It might be better suited to e.g. [Signal Processing](http://dsp.stackexchange.com/help/on-topic) or maybe [Cross-Validated](http://stats.stackexchange.com/help/on-topic) – jonrsharpe Jun 19 '14 at 11:16
  • Apologies, though I feel its the code implementation I'm getting wrong. I don't fully understand how the numpy functions are working. though if you feel its off topic I can of course move it – Leavenotrace Jun 19 '14 at 11:18
  • Is it that you think the implementation is incorrect? What are the results you were expecting? Have you tested the code - what works, what doesn't? A request to explain code would be off-topic here, and a general *"what's going wrong?"* is hard to answer - see [How to Ask](http://stackoverflow.com/help/how-to-ask) and consider providing a [minimal example](http://stackoverflow.com/help/mcve) – jonrsharpe Jun 19 '14 at 11:21
  • You're right, I'm new to stackoverflow and still getting used to it. I'm not asking anyone to explain how the numpy library works, I'm simpy saying I think tat might be where my error is and someone with a bit more experience might pick it up. I tried to address at the bottom that I have tested the code using a sin function and produced a peak corresponding to a single frequency. The white noise should give a flat power spectral denisty. basically it should look like the first half of the graph I have provided, leading me to believe I'm mirroring the data somewhere. – Leavenotrace Jun 19 '14 at 11:26
  • Your not getting some sort of boundary effect? If you signal and the thing your correlating don't completely overlap you will have smaller values at the edges as there are fewer values contributing. – Salix alba Jun 19 '14 at 13:12

1 Answers1

0

First,your function is written incorrect, you take fft from autocorr which is not your initial signal array. Funny thing, because of the nature of rounding error you will get whit noise like psd from it as well. Second, you calculating frequencies axis wrong, as they should be extended to N/(t*2) (which is Nyquist fr.). Instead, since your freq_axis array length is N, you retrieve N elements from you sig array, and because of that you just read negative frequencies with the same psd, which causes you the symmetry (taking log convert you variable to real and all the Fourier coefficients for negative freqs are just conjugates for positive ones). replacing your code with the following one give a perfectly good result:

sig = white(N,1,N/t)
(siglog,freq_axis)=ml.psd(sig,N,(N/t), detrend=ml.detrend_none,
   window=ml.window_hanning, noverlap=0, pad_to=None,
    sides='default', scale_by_freq=None)
plt.plot(freq_axis,np.log10(siglog))
plt.show()

matplotlib.mlab.psd result

don't forget to import

import matplotlib.mlab as ml
Andrey Debolskiy
  • 1