Engee documentation
Notebook

Construction of frequency response from experimental data

In this example, we will show how to plot the frequency response of an object from experimental data, either from a single output signal or from two signals - input and output.

Loading data

We have a file with one time signal.

In [ ]: 
Pkg.add(["Statistics", "DSP", "ControlSystems", "CSV"])
In [ ]: 
using CSV, DataFrames
d = CSV.read( "$(@__DIR__)/scope_data.csv", DataFrame );
plot( d.Time, d.Measurement, size=(500,300) )
Out[0]:

Let's take the average difference between time steps as the sampling frequency.

In [ ]: 
using Statistics
Fs = 1 / mean(diff(d.Time));    # Частота дискретизации
Fn = Fs/2;                      # Частота Найквиста (среза)
L = length( d.Time );           # Количество отсчетов в сигнале

Frequency characteristic of the signal

To translate a signal into the frequency domain, all you need is the function fft.

In [ ]: 
using FFTW
TF = fft( d.Measurement ) ./ L

Fl = round.(Int, round(L/2)+1)           # Количество точек в дискретном преобразовании Фурье
Fv = range(1/Fl, 1, length=Fl) .* Fn;    # Вектор частот (исключая первую точку – 0)
Iv = 1:length(Fv);                       # Вектор индексов

Let's also look at the smoothed spectrum - create an elliptic filter:

In [ ]: 
using DSP

Wp = (22, 30) ./ Fn; # Границы полосы пропускания
Ws = (20, 35) ./ Fn; # Границы полосы поглощения
Rp, Rs = 1, 30; # Пики пропускания и поглощения

n,Ws = ellipord( Wp, Ws, Rp, Rs);            # Найдем порядок фильтра с желаемыми характеристиками
designmethod = Elliptic( n, Rp, Rs );        # Создадим эллиптический фильтр
responsetype = Bandpass( Rp, Rs; fs=Fn );    # И полосной фильтр

# Сглаженный спектр
TF = fft( d.Measurement ) ./ L
TF1 = filt( digitalfilter(responsetype, designmethod), TF );

plot(
    plot( Fv, 20 .* log10.(abs.([TF[Iv] TF1[Iv]])), xscale=:log10, lw=[1 2] ),
    plot( Fv, angle.([TF[Iv] TF1[Iv]]) .* 180/pi, xscale=:log10, lw=[1 2] ),
    layout=(2,1)
)

plot!( ylabel=["Амплитуда (дБ)" "Фаза (град)"], leg=:false )
plot!( title=["ЛАФЧХ системы" ""], xlabel=["" "Частота (Гц)"] )
Out[0]:

The filter should be chosen in such a way that it is appropriate for the problem to be solved. Here we have implemented the filter without any attention to the result.

Frequency response of the filter

Now we are given an input signal and an output signal. We need to see the frequency response of the link that translates one signal into the other.

Our output signal is the same input signal passed through a filter with a transfer function $\frac{s+1}{s^2 + 0.8s + 1}$.

In [ ]: 
using ControlSystems
G = tf( [1,1], [1, 0.8, 1], Fs);
output, _, _ = ControlSystems.lsim( G, input );

Let's obtain the frequency response of the link that returns the output signal output with the input signal input:

In [ ]: 
FTinp = fft(input) ./ L;
FTout = fft(output) ./ L;
TF = FTout ./ FTinp;

# Добавим еще один график со сглаженным спектром
TF1 = filt( digitalfilter(responsetype, designmethod), TF );

plot(
    plot( Fv, 20 .* log10.(abs.([TF[Iv] TF1[Iv]])), title="Amplitude", ylabel="dB", xscale=:log10, lw=[1 2]),
    plot( Fv, angle.([TF[Iv] TF1[Iv]]) .* 180/pi, title="Phase", ylabel= "°", xscale=:log10, lw=[1 2]),
    layout=(2,1)
)

plot!( ylabel=["Амплитуда (дБ)" "Фаза (град)"], leg=:false )
plot!( title=["ЛАФЧХ системы" ""], xlabel=["" "Частота (Гц)"] )
Out[0]:

Conclusion

To obtain the frequency response of a process for which we are given the input and output signals, it is enough to divide one signal by the other after translating them into the frequency domain.