Engee documentation
Notebook

Construction of the frequency response based on experimental data

In this example, we will show how to construct the frequency response of an object based on experimental data, both from one output signal and two input and output signals.

Uploading 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 the time steps as the sampling frequency.

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

Frequency response of the signal

To translate the signal into the frequency domain, the function is sufficient 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 – let's create an elliptical 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 must be selected in such a way that it corresponds to the task being solved. Here we have implemented a filter without any attention to the result.

Frequency response of the filter

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

Our output signal will be the same input signal passed through a filter
with a transfer function. .

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

We obtain the frequency response of the link, which is at the input signal input returns the output signal output:

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 the process for which we are given the input and output signals, it is enough to divide one signal into another after their translation into the frequency domain.