恩吉的过滤功能
在本示例中,我们将研究指定信号的可能性、过滤信号的方法和改变信号频率的方法,并分析在与信号交互的各个阶段对信号功率谱密度的估计。
连接库
Pkg.add(["Noise", "DSP", "DigitalComm"])
using Plots; # Билиотека отрисовки графиков
using Noise; # Библиотека добавления шумов
using DigitalComm; # Библиотека систем связи
using DSP; # Библиотека цифровой обработки сигналов
using FFTW; # Библиотека преобразований Фурье
plotlyjs(); # Подключения метода отрисовки
设置正弦波
示例中的所有变换都将在该正弦波上进行。
fs = 5000; # Шагов в секунду
t = [0:1/fs:1;]; # Время
y = sin.(2*pi*10*t); # Задание синусоиды
plot(t,y)
让我们通过 periodogram
来估算信号的功率谱密度--这是一个周期图构造函数,基于对数据序列的傅里叶变换模的平方的计算。
物理学和信号处理中的功率谱密度(PSD)是描述信号功率随频率分布的函数,即单位频率间隔内的功率。
y_p = DSP.periodogram(y, onesided=true, nfft=length(y), fs=fs);
plot(freq(y_p)[1:20], power(y_p)[1:20], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
量化
信号转换的下一阶段是量化--将信号的参考值范围划分为有限数量的电平,并将这些值四舍五入到最接近的两个电平之一。
y_quant = quantization(y, 80, minv=-1, maxv=1); # Квантование исходной синусоиды
# Вывод результата
plot(y[1:250], label="исходный сигнал")
plot!(y_quant[1:250], label="квантованный сигнал по 80 уровням")
y_p = DSP.periodogram(y_quant, onesided=true, nfft=length(y_quant), fs=fs);
plot(freq(y_p)[1:20], power(y_p)[1:20], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
我们可以看到,这种转换不会影响信号的功率谱密度,也不会移动信号的频率。
数据减薄
我们接下来要对原始正弦波进行的操作是数据减薄。我们将从原始信号中每隔一秒取样一次。
y_d = y_quant[1:2:end]
t_new = t[1:2:end]
plot(t_new,y_d)
y_p = DSP.periodogram(y_d, onesided=true, nfft=length(y_d), fs=fs);
plot(freq(y_p)[1:20], power(y_p)[1:20], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
周期图显示,数据减薄会增加信号的频率。这是由于采样数量减少了两倍,振荡频率也相应增加了两倍。
噪声叠加
噪音*是阻碍我们从有用信号中获取信息的无用现象。
噪音具有随机性。
在数字信号处理中,**加性白高斯噪声(ABGN)**最为常用。
其特点如下
- 均匀的功率谱密度(因此称为白噪声);
- 时间值的正态分布(因此称为高斯);
- 与有用信号相加(因此称为加法);
- ABGS 在统计上独立于有用信号。
signal = addNoise(y_d, 20); # Добавление шумов
plot(t_new, signal, label="синусоида с шумом")
plot!(t_new, y_d, label="синусоида без шума")
值得注意的是,函数addNoise
定义了一个元组,我们只需要从中提取第一个向量。
下面的代码显示了输入addNoise
的数据类型、其输出以及从元组中提取第一个向量的结果。
typeof(y_d)
typeof(signal)
signal=signal[1];
typeof(signal)
Engee 中的信号过滤
在本例中,我们将介绍三种滤波器,但您可以在Engee文档中找到所需的任何滤波器。
我们将考虑的第一种方法是通过滤波器系数来指定滤波器,其形式为两个矢量a
和b
,然后通过传递函数的分子和分母的系数来实现滤波器的表示。
b
和分母a
的系数来实现滤波器的表示,使用函数库DSP
中的函数PolynomialRatio
。
注意
b
和a
可以指定为对象Polynomial
, 或
矢量,从高到低排列。
b = [0.75, 0.5];
a = [1, 0.2, 0.1];
f = PolynomialRatio(b,a);
Out_sig = filt(f,signal);
这个滤波器不是一个抽取滤波器。这可以通过比较滤波器输入和输出的大小来验证。
size(Out_sig,1)==size(signal,1)
plot(signal, label="исходный сигнал")
plot!(Out_sig, label="отфильтрованный сигнал")
指定相同滤波器的第二种方法是使用函数impresp
得出脉冲响应,然后使用函数conv
计算滤波器输出。
h = impresp(f);
Out_sig2 = conv(signal,h);
plot(signal, label="исходный сигнал")
plot!(Out_sig2, label="отфильтрованный сигнал")
让我们比较一下函数 filt
的输出和函数 conv
的输出。
plot(Out_sig[1:2500]-Out_sig2[1:2500])
从图中我们可以看到,函数filt
的输出和函数conv
的输出几乎完全重合。
我们要考虑的第三种滤波方法是设置窗口 FIR 滤波器。它与前几种滤波器不同,因为它相对于输入信号的延迟等于窗口大小。
responsetype = Lowpass(5; fs=30) # Определение полосы пропускания
designmethod = FIRWindow(hanning(128)) # Определение размеров окна
Out_fir = filt(digitalfilter(responsetype, designmethod), Out_sig)
plot(Out_sig, label="исходный сигнал")
plot!(Out_fir, label="исходный сигнал")
为了比较滤波器的输出,我们将使用韦尔奇周期图。
n=div(length(Out_sig),8);
y_out = welch_pgram(Out_sig, n, div(n,2), onesided=true, nfft=length(Out_sig), fs=fs);
n=div(length(Out_sig2),8);
y_in = welch_pgram(Out_sig2, n, div(n,2), onesided=true, nfft=length(Out_sig2), fs=fs);
n=div(length(Out_fir),8);
y_fir = welch_pgram(Out_fir, n, div(n,2), onesided=true, nfft=length(Out_fir), fs=fs);
plot(freq(y_in)[1:25], power(y_in)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", label="filt")
plot!(freq(y_out)[1:25], power(y_out)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", label="conv")
plot!(freq(y_fir)[1:25], power(y_fir)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", label="fir")
我们可以看到,滤波器输出的主要区别在于其功率谱密度。
改变信号频率
为了增加或减少采样次数,Engee 使用resample
命令,可以对信号执行upsampling
和downsampling
。
注:upsampling
和downsampling
通过关系式n//k
来指定,其中:
n
:要创建的点数;
k
每个点的近邻数。
# Downsampling
rs = resample(Out_fir, 1//2)
n=div(length(rs),8);
y_out = welch_pgram(rs, n, div(n,2), onesided=true, nfft=length(rs), fs=fs);
plot(freq(y_out)[1:25], power(y_out)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
我们可以看到,滤波信号的频率增加了一倍。
# Upsampling
rs2 = resample(Out_fir, 2//1)
n=div(length(rs2),8);
y_out = welch_pgram(rs2, n, div(n,2), onesided=true, nfft=length(rs2), fs=fs);
plot(freq(y_out)[1:25], power(y_out)[1:25], xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
我们可以看到,滤波信号的频率降低了一半。
println("Размер входа: ",size(Out_fir,1))
println("Размер при downsampling: ",size(rs,1))
println("Размер при upsampling: ",size(rs2,1))
plot(rs)
plot!(rs2)
在比较信号及其大小时也可以观察到同样的效果。在downsampling
时,信号的样本数会减少,而在upsampling
时,相对于原始信号,样本数反而会增加。
结论
在本例中,我们展示了与信号交互的可能性,并证明Engee完全适用于数字信号处理。