# Find the dominant low frequency

I have this signal, with an arbitrary x-axis, and I want to find the (visually) dominant frequency of the signal. There may be a high frequency component, but I'm interested in the component with a period of around 120 pixels.

First I implemented the Pisarenko Harmonic Decomposition, but only got a very high frequency out (0.04 Hz) (one bin, one second). The I tried with the Fourier transform, and when I plotted the frequency amplitude spectrum I also only got a range of very high frequencies.

Can someone help me solve this problem? What am I doing wrong? It feels like such a simple problem, but for some reason I just can't solve it.

The raw dataset is located here: http://pastebin.com/5ghKbBRm

• Is there a particular reason why you chose Pisarenko's method? It's kind of special-purpose, isn't it. I mean, I preach MUSIC/ESPRIT all over, but maybe a simple argmax(autocorrelation estimate) would do here? – Marcus Müller Jan 17 '17 at 10:21
• Also, I don't see your "very high frequency content"; could you elaborate? – Marcus Müller Jan 17 '17 at 10:22
• @MarcusMüller I first picked PHD just because it looked like it was doing what I needed. Maybe I was wrong. I'm also not sure there is a high frequency, but I just assumed there might be since that what's both my FFT and PHD returned – Markus Jan 17 '17 at 10:24
• I really don't see where your |FFT| shows high frequencies? Could you actually give me the coordinate in x/y axis values that you mean? – Marcus Müller Jan 17 '17 at 10:25
• @MarcusMüller Erm - I'm sorry. I must have been looking at the too long. I obviously talking complete nonsense. There are no high frequencies in my |FFT|, only low. As far as I can see they do not actually correspond to the frequency I'm looking for, but I think I have to re-evaluate my question. Sorry. However, could you tell me a little about MUSIC/ESPRIT and why I should use it? – Markus Jan 17 '17 at 10:29

Have a look at the Python code below. It calculates the auto-correlation of the data, after the DC part of the signal was removed. The autocorrelation shows nice peaks with much less noise than the original signal. Then, we detect the local maxima and measure the distance between them. This distance is then the estimated period of your signal:

import csv
import requests
import numpy as np
from io import BytesIO
from scipy import signal

url = 'http://pastebin.com/raw/5ghKbBRm'
response = requests.get(url)
A = np.genfromtxt(BytesIO(response.content)) # read the CSV file

# begin signal processing

A = A - np.mean(A)  # remove the DC part of the signal
plt.subplot(121)
plt.plot(A)

plt.subplot(122)
C = np.correlate(A, A, 'same')  # calculate autocorrelation to smooth out the noise
plt.plot(C)

plt.subplot(122)
P = signal.argrelextrema(C, np.greater)[0]  # caluclate local maxima (this step can be improved)
P = P[C[P]>15000] # filter maxima, where the autocorrelation values is too low
plt.stem(P, 30000*np.ones_like(P))  # plot where the detected peaks are

print ("Peak distances: ", np.diff(P))  # calculate distance between peaks and the mean of the distances
print ("Estimated period: ", np.mean(np.diff(P)))


Program output:

Peak distances:  [117 118 117 117 117 117 118 117]
Estimated period:  117.25


There is room for improvement. For example, the detection of local maxima is not very robust and the threshold was set manually, but this should enable you to find out the rest on your own.