Frequency Estimation

In [1]:

    

from matplotlib import pyplot
from bokeh import plotting as plt
from bokeh.io import output_notebook
from scipy import interpolate
from scipy.signal import argrelextrema
from scipy.signal import savgol_filter
from scipy.stats import signaltonoise

from IPython.display import display, HTML

import os
import pandas as pd
import numpy as np
import numba as nb
import sys


    

In [2]:

    

%matplotlib inline  
output_notebook()


    

    






    

        
        Loading BokehJS ...
    






    

In [3]:

    

%%javascript
IPython.OutputArea.auto_scroll_threshold = 9999;


    

In [4]:

    

# constants
π = np.pi


    

In [5]:

    

def gen_sine_signal(
    t: np.array, A: float, f: float, φ: float
) -> pd.Series:
    """
    φ input in degree unit

    :param t:
    :type t:
    :param A:
    :type A:
    :param f:
    :type f:
    :param φ:
    :type φ:
    :return:
    :rtype:
    """
    ω = 2*np.pi*f
    return pd.Series(A*np.sin(ω*t + np.deg2rad(φ)), index=t)


    

In [6]:

    

def calculate_frequency(
    data: np.array, fs: float, ns:float, low_pass: float=None
):
    fs = np.round(fs).astype(int)
    time_interval = 1/fs
    ts = ns//fs
    limit = int(ts*fs)
    
    data_freq = np.fft.fft(data[:limit].copy())
    N = data_freq.size
    
    # gives us a list of frequencies for the FFT
    w = np.fft.fftfreq(N, time_interval)
    
    if low_pass:
        ipos = np.where(np.all((w > 0, w < low_pass), axis=0))
    else:
        ipos = np.where(w > 0)
        
    freqs = w[ipos]  # only look at positive frequencies
    mags = abs(data_freq[ipos])  # magnitude spectrum
    
    # _freq = freqs[mags.tolist().index(max(mags))]
    return freqs[np.where(mags == max(mags))[0]][0]


    

In [7]:

    

def plot_data(df: pd.DataFrame, fs: int):
    # charts of filtered data
    for i in range(len(df.keys())):
        # data
        se = df.loc[:, i].reset_index(drop=True)

        # plot
        pf = plt.figure(
            width=800, height=200, 
            x_axis_label="TIME (s)",
            y_axis_label="Force (N)",
            title='Signal Data - CH ' + str(i)
        )
        pf.line(
            se.index, 
            se.values, 
            legend='cycles',
        )

        pf.line(
            se.index[:fs], 
            se.values[:fs], 
            legend='cycle with limit (%s)' % fs,
            color='red'
        )

        plt.show(pf)


    

In [8]:

    

# GENERAL INFO
channels = 1
φs = [45]  # phase shift of each channel

f = 25

ns = 680
fs = 680

t = np.linspace(0, ns//fs, ns)

df = pd.DataFrame(np.zeros((ns, channels)), index=t)

for i in range(channels):
    φ = φs[i]
    title = 'Phase of channel %s: %s (deg)' % (i, φ)
    df[i] = gen_sine_signal(t, A=1, f=f, φ=φ)


    

In [9]:

    

# plot the data channel and data with between 0 until "fs" value
plot_data(df, fs)


    

In [10]:

    

_text = '(CH %s | fs %s | ns %s): %s Hz'

display(HTML(
    '<h4>Frequency estimation - Pure signal - %s Hz</h4>' % f
))

for k, v in df.items():
    fq = calculate_frequency(
        v.values, fs=fs, ns=ns
    )
    print(_text % (k, fs, ns, fq))


    

    





Frequency estimation - Pure signal - 25 Hz






    




(CH 0 | fs 680 | ns 680): 25.0 Hz

In [11]:

    

df_noise_data = df.copy()

_text = '(CH %s | fs %s | ns %s): %s Hz'

display(HTML('''
    <h4>Frequency estimation - 
    Noisy signal (low amplitude noise) -
    %s Hz</h4>
    ''' % f
))

for k, v in df.items():
    df_noise_data[k] = v+gen_sine_signal(t, A=1, f=900, φ=0)
    
    fq = calculate_frequency(
        df_noise_data[k], fs=fs, ns=ns
    )
    print(_text % (k, fs, ns, fq))

plot_data(df_noise_data, fs)


    

    






    
Frequency estimation - 
    Noisy signal (low amplitude noise) -
    25 Hz
    






    




(CH 0 | fs 680 | ns 680): 25.0 Hz






    







    

In [12]:

    

df_noise_data = df.copy()

_text = '(CH %s | fs %s | ns %s): %s Hz'

display(HTML('''
    <h4>Frequency estimation - 
    Noisy signal (high amplitude noise):
    %s Hz</h4>
    ''' % f
))

# apply noise
for k, v in df.items():
    df_noise_data[k] = v+gen_sine_signal(t, A=2, f=200, φ=0)

# plot data with noise
plot_data(df_noise_data, fs)

# frequency estimation
for k, v in df.items():
    fq = calculate_frequency(
        df_noise_data[k], fs=fs, ns=ns
    )
    print(_text % (k, fs, ns, fq))    


display(HTML('''
<h5>Estimation with frequency filter low pass</h5>
'''))
# frequency estimation using filter by low pass
for k, v in df.items():
    fq = calculate_frequency(
        df_noise_data[k], fs=fs, ns=ns, low_pass=60
    )
    print(_text % (k, fs, ns, fq))


    

    






    
Frequency estimation - 
    Noisy signal (high amplitude noise):
    25 Hz
    






    







    








    




(CH 0 | fs 680 | ns 680): 200.0 Hz






    






Estimation with frequency filter low pass







    




(CH 0 | fs 680 | ns 680): 25.0 Hz