Find the Average Frequency

33 Views Asked by At

The average frequency of a function $x(t)$ is defined as:

$$ \langle \omega \rangle = \frac{ \int_{-\infty}^{\infty} \omega \cdot | X(\omega) |^{2} d \omega }{\int_{-\infty}^{\infty} | X(\omega) |^{2} d \omega} $$

where $X(\omega)$ is the Fourier transform of $x(t)$. I want to avoid taking the Fourier transform, so instead, according to Leon Cohen in his book Time Frequency Analysis, I can compute the average frequency via:

$$ \langle \omega \rangle = \frac{1}{j} \cdot \frac{\int_{-\infty}^{\infty} x^{\ast}(t) \cdot \frac{\partial }{\partial t} x(t) dt}{ \int_{-\infty}^{\infty} |x(t)|^{2} dt }$$

When $x(t) = A(t) e^{j \phi(t)}$, this reduces to:

$$ \langle \omega \rangle = \frac{ \int_{-\infty}^{\infty} \frac{d}{dt} \phi(t) \cdot A^{2}(t) dt }{ \int_{-\infty}^{\infty} |x(t)|^{2} dt }$$

Now I am trying to implement this in Python, but I am having no luck. I do not know what I am doing wrong, and wondering if anyone has any idea. My code is below -- (pieces stolen from around the internet):

import numpy
import matplotlib.pyplot as plt
import scipy.integrate as integrate
import scipy.signal as signal
import quadpy

def cosine_wave(f,fs,phase,num_cycles):
    t = numpy.arange(0,num_cycles*1/f-1/fs,1/fs) 
    g = numpy.cos(2*numpy.pi*f*t+phase) 
    return t,g 


def derivative(f,a,method='central',h=0.01):
    if method == 'central':
        return (f(a + h) - f(a - h))/(2*h)
    elif method == 'forward':
        return (f(a + h) - f(a))/h
    elif method == 'backward':
        return (f(a) - f(a - h))/h
    else:
        raise ValueError("Method must be 'central', 'forward' or 'backward'.")
    return

class signal_moment_calculator:

    def __init__(self, signal, time_samples, sampling_rate):
        self.signal = numpy.copy(signal)
        self.time = numpy.copy(time_samples)
        self.fs = numpy.copy(sampling_rate)
        self.signal_energy = integrate.simps(y=abs(self.signal)**2.0, x=self.time, dx=abs(self.time[1]-self.time[0]))
        self.signal_derivative = derivative(self.sample, time_samples, method='central', h = abs(time_samples[1] - time_samples[0]))

        self.signal_amplitude = abs(self.signal)
        self.signal_angles = numpy.unwrap(numpy.angle(self.signal)) 
        self.signal_angle_derivative = derivative(self.sample_angle, time_samples, method='central', h = abs(time_samples[1] - time_samples[0]))

    def sample(self, t_val):
        idx = []
        for t_v in t_val:
            idx.append(numpy.abs(self.time - t_v).argmin())     
        return numpy.asarray((self.signal[idx]))

    def sample_angle(self, t_val):
        idx = []
        for t_v in t_val:
            idx.append(numpy.abs(self.time - t_v).argmin())     
        return numpy.asarray((self.signal_angles[idx]))

    def mean_freq_calculator(self, t_val):
        idx = []
        for t_v in t_val:
            idx.append(numpy.abs(self.time - t_v).argmin())     
        return (1.0/1j) * numpy.conjugate(self.signal[idx]) * self.signal_derivative[idx]
        #return self.signal_angle_derivative[idx] * (abs(self.signal_amplitude[idx])**2.0)
        #return self.signal_angle_derivative[idx] * 2.0*(abs(self.signal_real[idx])**2.0)

    def calculate_freq_moments(self):
        mean,error = quadpy.quad(self.mean_freq_calculator, self.time[0], self.time[-1], epsabs=1e-4, epsrel=1e-4, limit=5000)
        mean /= self.signal_energy
        self.average_freq_ = mean
        return

#==============TEST===============#
freq = 1000
fs = 20*freq
num_cycles = 10
t_cos, x_cos = cosine_wave(freq,fs,0.0,num_cycles)
x_analytic = signal.hilbert(x_cos)

moment_calculator = signal_moment_calculator(x_analytic, t_cos, fs)
moment_calculator.calculate_freq_moments()
print("True Average Frequency: ", freq, "Average Frequency Calculated: ", moment_calculator.average_freq)    

No matter what the frequency is, the sampling rate, or the number of cycles, or which of the three methods (that are commented out), the numbers are nowhere near each other.

Questions:

1. What am I doing wrong? 

2. Is there another (easier) way to calculate the average frequency without resorting to taking the Fourier transform?

Thanks!