根据实验数据构建频率响应

在本示例中，我们将展示如何根据实验数据绘制物体的频率响应，实验数据可以是单一输出信号，也可以是两个信号（输入和输出）。

加载数据

我们有一个包含一个时间信号的文件。

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) )

我们将各时间步之间的平均差作为采样频率。

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

信号的频率特性

要将信号转换到频域，只需使用函数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);                       # Вектор индексов

我们也来看看平滑频谱--创建一个椭圆滤波器：

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=["" "Частота (Гц)"] )

滤波器的选择应与要解决的问题相匹配。在这里，我们在不考虑结果的情况下使用了滤波器。

滤波器的频率响应

现在我们得到一个输入信号和一个输出信号。我们需要查看将一个信号转换成另一个信号的链路的频率响应。

我们的输出信号是相同的输入信号通过一个具有传递函数的滤波器
的滤波器，其传递函数为。

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

让我们来获取输出信号output 与输入信号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=["" "Частота (Гц)"] )

结论

要获得给定输入和输出信号的过程的频率响应，只需将一个信号除以另一个信号，然后将它们转换到频域即可。