Документация Engee
Notebook

Построение частотной характеристики по экспериментальным данным

В этом примере мы покажем, как построить частотную характеристику объекта по экспериментальным данным – как по одному выходному сигналу, так и по двум – входному и выходному.

Загрузка данных

У нас есть файл с одним временным сигналом.

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]:

В качестве частоты дискретизации возьмем усредненную разницу между временными шагами.

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

Частотная характеристика сигнала

Чтобы перевести сигнал в частотную область, достаточно функции 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);                       # Вектор индексов

Посмотрим также на сглаженный спектр – создадим эллиптический фильтр:

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]:

Фильтр должен подобрать таким образом, чтобы он соответствовал решаемой задаче. Здесь мы реализовали фильтр без какого-либо внимания к результату.

Частотная характеристика фильтра

Теперь нам даны входной и выходной сигнал. Нужно увидеть частотную характеристику звена, которое переводит один сигнал в другой.

Нашим выходным сигналом будет тот же входной сигнал, пропущенный через фильтр
с передаточной функцией .

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

Получим частотную характеристику звена, которое при входном сигнале input возвращает выходной сигнал 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]:

Заключение

Чтобы получить частотную характеристику процесса, для которого нам даны входной и выходной сигналы, достаточно разделить один сигнал на другой после их перевода в частотную область.