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));    # Sampling rate
Fn = Fs/2;                      # Nyquist frequency (cutoff)
L = length( d.Time );           # The number of samples in the signal

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)           # The number of points in the discrete Fourier transform
Fv = range(1/Fl, 1, length=Fl) .* Fn;    # Frequency vector (excluding the first point – 0)
Iv = 1:length(Fv);                       # Index vector

Let's also look at the smoothed spectrum – let's create an elliptical filter.:

In [ ]: 
using DSP

Wp = (22, 30) ./ Fn; # Bandwidth limits
Ws = (20, 35) ./ Fn; # Absorption band boundaries
Rp, Rs = 1, 30; # Transmission and absorption peaks

n,Ws = ellipord( Wp, Ws, Rp, Rs);            # Let's find the filter order with the desired characteristics
designmethod = Elliptic( n, Rp, Rs );        # Creating an elliptical filter
responsetype = Bandpass( Rp, Rs; fs=Fn );    # And a bandpass filter

# Smoothed spectrum
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=["Amplitude (dB)" "Phase (hail)"], leg=:false )
plot!( title=["LFCH systems" ""], xlabel=["" "Frequency (Hz)"] )
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;

# Let's add another graph with a smoothed spectrum.
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=["Amplitude (dB)" "Phase (hail)"], leg=:false )
plot!( title=["LFCH systems" ""], xlabel=["" "Frequency (Hz)"] )
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.