基于实验数据的频率响应构造

在这个例子中，我们将展示如何根据一个输出信号和两个输入和输出信号的实验数据构造对象的频率响应。

上传数据

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

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=["振幅(dB)" "阶段(冰雹)"], leg=:false )
plot!( title=["LFCH系统" ""], xlabel=["" "频率(Hz)"] )

筛选器的选择必须与正在解决的任务相对应。在这里，我们已经实现了一个过滤器，没有任何关注结果。

滤波器的频率响应

现在给出了一个输入和一个输出信号。您需要查看将一个信号转换为另一个信号的链路的频率响应。

我们的输出信号将是通过滤波器的相同输入信号
具有传递函数。 .

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

我们得到链路的频率响应，它在输入信号处 input 返回输出信号 output:

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=["振幅(dB)" "阶段(冰雹)"], leg=:false )
plot!( title=["LFCH系统" ""], xlabel=["" "频率(Hz)"] )

结论

为了获得给出输入和输出信号的过程的频率响应，将一个信号转换到频域后分成另一个信号就足够了。