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.
Pkg.add(["Statistics", "DSP", "ControlSystems", "CSV"])
using CSV, DataFrames
d = CSV.read( "$(@__DIR__)/scope_data.csv", DataFrame );
plot( d.Time, d.Measurement, size=(500,300) )
Let's take the average difference between time steps as the sampling frequency.
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
.
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:
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=["" "Частота (Гц)"] )
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 .
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
:
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=["" "Частота (Гц)"] )
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.