CODEWITHSUNDEEP
frequencycorrelationspectrumpearson
branny12000
branny12000

Reputation: 185

comparing two frequency spectra

I'm trying to compare two frequency spectra but I am confused over a number of points.

One device samples at 40 Hz the other at 100 Hz and so I'm not sure if I need to take this into account. Anyway I have produced frequency spectra from both devices and now I wish to compare these. How can I do correlation at each point so that I get pearson correlations at each point. I know how to do an overall one of course but I want to see points of strong correlation and those less strong?

Upvotes: 0

Views: 1576

Answers (1)

Mitch
Mitch

Reputation: 81

If you are calculating power spectral densities P(f), then it doesn't matter how your original signal x(t) is sampled. You can directly and quantitatuively compare both spectra. To make sure that you have calculated the spectral densities you can explicitly check Parsevals theorem:

$ \int P(f) df = \int x(t)^2 dt $

Of course you have to think about which frequencies are actually evaluated Remember that a fft gives you frequencies from f = 1/T until or below the Nyquist frequency f_ny = 1/(2 dt) depending on the number of samples in x(t) being even or odd.

Here's a python example code for psd

def psd(x,dt=1.):

"""Computes one-sided power spectral density of x.

PSD estimated via abs**2 of Fourier transform of x
Takes care of even or odd number of elements in x:
- if x is even both f=0 and Nyquist freq. appear once
- if x is odd f=0 appears once and Nyquist freq. does not appear

Note that there are no tapers applied: This may lead to leakage!

Parseval's theorem (Variance of time series equal to integral over PSD) holds and can be checked via
print ( np.var(x), sum(Px*f[1]) )
Accordingly, the etsimated PSD is independent of time series length

Author/date: M. von Papen / 16.03.2017
"""

N = np.size(x)
xf = np.fft.fft(x)
Px = abs(xf)**2./N*dt
f  = np.arange(N/2+1)/(N*dt)

if np.mod(N,2) == 0:
    Px[1:N/2] = 2.*Px[1:N/2]
else:
    Px[1:N/2+1] = 2.*Px[1:N/2+1]

# Take one-sided spectrum
Px = Px[0:N/2+1]

return Px, f

Upvotes: 1

Related Questions