Syslab使用手册
⋮

# cwtfilterbank

连续小波变换滤波器组

函数库: TyWavelet

# 描述

使用 cwtfilterbank 创建连续小波变换 (CWT) 滤波器组。滤波器组中使用的默认小波是解析 Morse(3, 60) 小波。您可以改变 Morse 小波的时间带宽和对称参数，以根据您的需要调整 Morse 小波。您也可以使用分析 Morlet(Gabor) 小波或凹凸小波。在时频分析多个信号时，为了提高计算效率，可以预先计算滤波器一次，然后将滤波器组作为输入传递给 cwt。使用滤波器组，您可以在时间和频率上可视化小波。您还可以创建具有特定频率或周期范围的滤波器组，并测量 3dB 带宽。您可以确定滤波器组中小波的质量因子。

# 构造

fb = cwtfilterbank( )
fb = cwtfilterbank(; Name = Value)

# 说明

fb = cwtfilterbank( ) 创建连续小波变换 (CWT) 滤波器组 fb。滤波器被归一化，使得所有通带的峰值幅度近似等于 2。默认滤波器组是为具有 1024 个样本的信号设计的。默认滤波器组使用解析 Morse(3, 60) 小波。滤波器组使用默认比例：每个倍频程大约有 10 个小波带通滤波器(每个倍频程有 10 个声音)。最高频率通带被设计为使得幅度下降到奈奎斯特频率处的峰值的一半。

正如所实现的，CWT 使用 L1 归一化。利用 L1 归一化，不同尺度的等幅值振荡分量在 CWT 中具有相等的幅度。L1 归一化提供了信号的更精确的表示。振荡分量的幅值与相应的小波系数的幅值一致。

fb 可以用作 cwt 的输入。

fb = cwtfilterbank(; Name = Value) 使用一个或多个名称值参数创建具有 Properties 的 CWT 筛选器组 fb。属性可以按任何顺序指定为 Name1 = Value1，...，NameN = ValueN。

TIP

不能更改现有筛选器组的属性值。例如，如果您有一个 SignalLength 为 2000 的滤波器组 fb，则必须创建第二个滤波器组 fb2 来处理具有2001 个样本的信号。不能为 fb 指定不同的 SignalLength。

# 属性

SignalLength - 信号长度
1024 (默认) | 正整数 ≥ 4

信号的长度，指定为正整数。信号必须至少有四个样本。

示例： fb = cwtfilterbank(SignalLength = 1700)

Wavelet - 解析小波
"Morse" (默认值) | "amor" | "bump"

滤波器组中使用的解析小波，指定为 "Morse"、"amor" 或 "bump"。这些字符串分别指定分析 Morse、Morlet (Gabor) 和凹凸小波。默认的小波是解析 Morse(3, 60) 小波。

默认情况下，对于 Morse 小波，频率响应在奈奎斯特衰减到峰值幅度的 50%。对于 Morlet 和 bump 小波，频率响应衰减到峰值幅度的 10%。您可以通过设置滤波器组 FrequencyLimits 属性来更改衰减百分比。

对于 Morse 小波，也可以使用 TimeBandwidth 或 WaveletParameters 属性对小波进行参数化。

示例： fb = cwtfilterbank(SignalLength = 1700, Wavelet = "bump")

VoicesPerOctave - 每个倍频程的声音数
10 (默认值) | 介于 1 和 48 之间的整数

用于 CWT 的每个倍频程的声音数，指定为 1 到 48 之间的整数。CWT 音阶是使用指定的每八度音的数量来离散化的。小波在频率和时间上的能量分布自动决定了最小和最大尺度。

SamplingFrequency - 采样频率(赫兹)
1 (默认值) | 正标量

采样频率(赫兹)，指定为正标量。如果未指定，频率以 cycles/sample 为单位，奈奎斯特频率为 ½。要以句点指定小数位数，请使用 SamplingPeriod 和 PeriodLimits 名称值参数。

不能同时指定 SamplingFrequency 和 SamplingPeriod 属性。

示例： fb = cwtfilterbank(SamplingFrequency = 5,Wavelet = "amor")

FrequencyLimits - 频率限制
二元标量向量

小波滤波器组的频率极限，指定为具有正严格递增项的两元素向量。

  • 第一个元素指定最低峰值通带频率。频率必须大于或等于以赫兹为单位的小波峰值频率和两个时间标准差除以信号长度的乘积。
  • 第二个元素指定最高峰值通带频率。高频限值必须小于或等于奈奎斯特。
  • 高频极限 fMax 与低频极限 fMin 之比的以 2 为底的对数必须大于或等于 1/NV，其中 NV 是每八度音阶的声音数量：

log2(fMax/fMin)≥1/NV

如果指定的频率限制超出了允许的范围，cwtfilterbank 会将限制截断为最小值和最大值。

如果在滤波器组中使用采样周期，则不能指定 FrequencyLimits 属性

示例： 如果 fb = cwtfilterbank(SignalLength = 1000, SamplingFrequency = 1000, FrequencyLimits = [90 100])，则 log2(100/90) ≥ 1/fb.VoicesPerOctave。

SamplingPeriod - 采样周期
持续时间标量

采样周期，指定为持续时间标量。不能同时指定 SamplingFrequency 和 SamplingPeriod 属性。

示例： fb = cwtfilterbank(SamplingPeriod = cwtSec(0.5))

注意 cwtfilterbank 用到以时间为单位的持续时间时，可用 cwtSec-秒、cwtMin-分钟、cwtHr-小时

PeriodLimits-周期限制
双元素持续时间阵列

小波滤波器组的周期极限，指定为具有正严格递增项的两元素持续时间阵列。

  • PeriodLimits 的第一个元素指定最大峰值通带频率，并且必须大于或等于 SamplingPeriod 的两倍。
  • 最大周期不能超过信号长度除以小波的两个时间标准差和小波峰值频率的乘积。
  • 最小周期 minP 与最大周期 maxP 之比的以 2 为底的对数必须小于或等于-1/NV，其中 NV 是每八度音阶的声音数量：

log2(minP/maxP) ≤ -1/NV

如果指定的周期限制超出了允许的范围，cwtfilterbank 会将限制截断为最小值和最大值。

如果在滤波器组中使用采样频率，则不能指定 PeriodLimits 属性。

示例：如果 fb = cwtfilterbank(SignalLength=1000,SamplingPeriod= cwtSec(0.1),PeriodLimits=[cwtSec(0.2) cwtSec(3)])，则log2(0.2/3) ≤ -1/fb.VoicesPerOctave。 数据类型：持续时间

TimeBandwidth - Morse小波的时间带宽乘积
60 (默认值) | 正标量

Morse 小波的时间带宽乘积，指定为大于或等于 3 且小于或等于 120 的正标量。Morse 小波的对称性 (gamma) 固定为 3。此属性仅在 Wavelet 为 "Morse" 时有效。

时间带宽乘积越大，小波在时间上越分散，小波在频率上越窄。Morse 小波在时间上的标准偏差大约为 sqrt(TimeBandwidth/2)。频率的标准偏差约为 1/2 × sqrt(2/TimeBandwidth)。

不能同时指定 TimeBandwidth 和 WaveletParameters 属性。

在 Morse 小波的表示法中，时间带宽是 。

示例： TimeBandwidth = 20

WaveletParameters - Morse小波参数
[3, 60] (默认) | 标量的两元素向量

Morse 小波参数，指定为两元素向量。第一个元素是对称参数 (gamma)，该参数必须大于或等于 1。第二个元素是时间带宽乘积，它必须大于或等于 gamma。时间带宽乘积与伽玛的比值不能超过 40。

当 gamma 等于 3 时，Morse 小波在频域中是完全对称的。偏斜度等于 0。gamma 值大于 3 会导致正偏度，而 gamma 值小于 3 会导致负偏度。

WaveletParameters 和 TimeBandwidth 名称值参数不能同时指定。

示例： fb = cwtfilterbank(WaveletParameters = [4, 20])

Boundary - 边界延伸
"reflection" (默认) | "periodic"

信号的边界扩展，指定为 "reflection" 或 "periodic"。确定在边界处如何处理数据。

TIP

如果要使用双帧或近似合成滤波器反转CWT，请将 Boundary 设置为 "periodic"。

示例： fb = cwtfilterbank(Boundary="periodic")

# 对象函数

名称 说明
wt 带滤波器组的连续小波变换
freqz CWT 滤波器组频率响应
timeSpectrum 时间平均小波频谱
scaleSpectrum 尺度平均小波频谱
wavelets 连续小波变换滤波器组时域小波
scales CWT 滤波器组刻度
waveletsupport CWT 滤波器组时间支持
qfactor CWT 滤波器组品质因数
powerbw CWT 滤波器组 3 dB 带宽
centerFrequencies CWT滤波器组带通中心频率
centerPeriods CWT滤波器组带通中心周期

# 示例

连续小波变换滤波器组

创建一个连续的小波变换滤波器组。

using TyWavelet

fb = cwtfilterbank()

cwtfilterbank - 属性:
VoicesPerOctave: 10
Wavelet: Morse
SamplingFrequency: 1
SamplingPeriod: nothing
PeriodLimits: cwtTime[]
SignalLength: 1024
FrequencyLimits: Real[]
TimeBandwidth: 60
WaveletParameters: Real[]
Boundary: reflection

绘制幅度频率响应图。

freqz(fb; plotfig=true)
连续小波变换滤波器组的频率分辨率

创建两个频率分别为 16 和 64 Hz 的正弦波。数据采样频率为 1000 Hz。绘制信号图。

using TyWavelet
using TyPlot

Fs = 1e3
t = collect(0:(1 / Fs):(1 - 1 / Fs))
x =
    cos.(2 * pi * 64 .* t) .* (t .>= 0.1 .&& t .< 0.3) .+
    sin.(2 * pi * 16 .* t) .* (t .>= 0.5 .&& t .< 0.9);
plot(t, x)
title("Signal")
xlabel("Time (s)")
ylabel("Amplitude")

为信号创建一个 CWT 滤波器组。绘制滤波器组中各小波的频率响应图。

fb = cwtfilterbank(; SignalLength=length(t), SamplingFrequency=Fs)
figure()
freqz(fb; plotfig=true)
title("Frequency Responses — Morse (3,60) Wavelet")

解析 Morse(3,60) 小波是滤波器库中的默认小波。该小波的时间带宽乘积等于 60。创建与第一个滤波器组相同的第二个滤波器组，但使用解析 Morse(3,5) 小波。绘制第二个滤波器组中各小波的频率响应图。

fb3x5 = cwtfilterbank(; SignalLength=length(t), SamplingFrequency=Fs, TimeBandwidth=5)
figure()
freqz(fb3x5; plotfig=true)
title("Frequency Responses — Morse (3,5) Wavelet")

请注意，频率响应比第一个滤波器组更宽。Morse(3,60) 小波的频率定位比 Morse(3,5) 小波更好。将每个滤波器组应用于信号，并绘制得到的频谱图。观察发现，Morse(3,60) 小波比 Morse(3,5) 小波具有更好的频率分辨率。

figure()
cwt(x; FilterBank=fb, plotfig=true)
title("Magnitude Scalogram — Morse (3,60)")

figure()
cwt(x; FilterBank=fb3x5, plotfig=true)
title("Magnitude Scalogram — Morse (3,5)")
正弦波和小波系数幅值

本例说明信号中振荡成分的幅值与相应小波系数的幅值一致。

创建一个由两个在时间上不相连的正弦波组成的信号。其中一个正弦波的频率为 32 Hz，幅值等于 1；另一个正弦波的频率为 64 Hz，幅值等于 2。绘制信号图。

using TyWavelet
using TyPlot

frq1 = 32
amp1 = 1
frq2 = 64
amp2 = 2
Fs = 1e3
t = collect(0:(1 / Fs):1)
x =
    amp1 .* sin.(2 * pi * frq1 .* t) .* (t .>= 0.1 .&& t .< 0.3) +
    amp2 .* sin.(2 * pi * frq2 .* t) .* (t .> 0.6 .&& t .< 0.9)
plot(t, x)
grid("on")
xlabel("Time (sec)")
ylabel("Amplitude")
title("Signal")

创建可应用于信号的 CWT 滤波器组。由于信号分量的频率是已知的，因此将滤波器组的频率限制设置为包括已知频率的较窄范围。为了确认范围，请绘制滤波器组的幅频响应图。

fb = cwtfilterbank(; SignalLength=length(x), SamplingFrequency=Fs, FrequencyLimits=[20 100])
freqz(fb; plotfig=true)

使用 cwt 和滤波器组绘制信号的频谱图。

figure()
cwt(x; FilterBank=fb, plotfig=true)
广义 Morse 和解析 Morlet 小波

本例展示了如何通过改变广义 Morse 小波的时间带宽参数来逼近解析莫里特小波。

广义 Morse 小波是精确解析小波的一个系列。Morse 小波有两个参数，即对称性和时间带宽乘积。您可以通过改变这些参数来获得具有不同特性和行为的解析小波。

加载 1995 年神户地震期间记录的地震仪数据。这些数据是 1995 年 1 月 16 日在澳大利亚霍巴特塔斯马尼亚大学记录的地震仪（垂直加速度，纳米/平方秒）测量数据，从 20:56:51 （格林尼治标准时间）开始，以 1 秒为间隔持续 51 分钟。创建默认设置的 CWT 滤波器库，并将其应用于数据。使用滤波器组生成示踪图。

using TyWavelet
using TyPlot
using TyBase

pkg_dir = pkgdir(TyWavelet)
source_path = pkg_dir * "/examples/Resources/kobe.mat"
y = load(source_path)

fb = cwtfilterbank(; SignalLength=length(kobe), SamplingFrequency=1)
cwt(kobe; FilterBank=fb, plotfig=true)

在 10 mHz 至 100 mHz 的频率范围内，小波系数的幅度较大。创建一个新的滤波器组，将频率限制设置为这些值。生成频谱图。

fb2 = cwtfilterbank(;
    SignalLength=length(kobe), SamplingFrequency=1, FrequencyLimits=[1e-2 1e-1]
)
cwt(kobe; FilterBank=fb2, plotfig=true)
title("Default Morse (3,60) Wavelet")

默认情况下，cwtfilterbank 使用 Morse(3,60) 小波。使用具有相同频率限制的解析莫列特小波创建滤波器组。生成一个频谱图，并与 Morse(3,60) 小波生成的频谱图进行比较。

fbMorlet = cwtfilterbank(;
    SignalLength=length(kobe),
    SamplingFrequency=1,
    FrequencyLimits=[1e-2 1e-1],
    Wavelet="amor",
)
cwt(kobe; FilterBank=fbMorlet, plotfig=true)
title("Analytic Morlet Wavelet")

Morlet 小波的频率局部性不如 (3,60) Morse 小波。不过，通过改变时间带宽乘积，您可以创建与 Morlet 小波特性相似的 Morse 小波。

使用 Morse 小波创建滤波器组，时间带宽值为 30 [2]，频率限制如上。生成地震仪数据的频谱图。请注意，频率上的涂抹现象与 Morlet 结果几乎相同。

fbMorse = cwtfilterbank(;
    SignalLength=length(kobe),
    SamplingFrequency=1,
    FrequencyLimits=[1e-2 1e-1],
    TimeBandwidth=30,
)
cwt(kobe; FilterBank=fbMorse, plotfig=true)
title("Morse (3,30)")

现在检查与 fbMorlet 和 fbMorse 滤波器组相关的小波。从这两个滤波器组中获取小波中心频率、滤波器频率响应和时域小波。确认中心频率几乎相同。

cfMorlet = centerFrequencies(fbMorlet)
frMorlet, fMorlet = freqz(fbMorlet)
wvMorlet, tMorlet = wavelets(fbMorlet)
cfMorse = centerFrequencies(fbMorse)
frMorse, fMorse = freqz(fbMorse)
wvMorse, tMorse = wavelets(fbMorse)
println("Number of Center Frequencies: " * string(length(cfMorlet)))

Number of Center Frequencies: 34

println("Maximum difference: " * string(maximum(abs.(cfMorlet - cfMorse))))

Maximum difference: 1.3877787807814457e-17

每个滤波器组都包含相同数量的小波。选择一个中心频率，绘制每个滤波器组中相关滤波器的频率响应。确认响应几乎相同。

wv = 13
figure()
plot(fMorlet, frMorlet[wv, :])
hold("on")
plot(fMorse, frMorse[wv, :])
grid("on")
hold("off")
title("Frequency Response")
xlabel("Frequency")
ylabel("Amplitude")
legend("Morlet", "Morse (3,30)")

绘制与相同中心频率相关的时域小波。确认它们几乎完全相同。

figure()
subplot(2, 1, 1)
plot(tMorlet, real(wvMorlet[wv, :]))
hold("on")
plot(tMorse, real(wvMorse[wv, :]))
grid("on")
hold("off")
title("Real")
legend("Morlet", "Morse (3,30)")
xlim([-100 100])
subplot(2, 1, 2)
plot(tMorlet, imag(wvMorlet[wv, :]))
hold("on")
plot(tMorse, imag(wvMorse[wv, :]))
grid("on")
hold("off")
title("Imaginary")
legend("Morlet", "Morse (3,30)")
xlim([-100 100])
更改时间 - 带宽积

这个例子说明，增加时间-带宽乘积 Morse 小波的包络线下会产生更多振荡。增加 会缩小小波的频率。

创建两个滤波器组。一个滤波器组的默认 TimeBandwidth 值为 60。第二个滤波器组的 TimeBandwidth 值为 10。两个滤波器组的 SignalLength 都是 4096 个样本。

using TyWavelet
using TyPlot
using TyMath

sigLen = 4096
fb60 = cwtfilterbank(; SignalLength=sigLen)
fb10 = cwtfilterbank(; SignalLength=sigLen, TimeBandwidth=10)

获取滤波器组的时域小波。

psi60, t = wavelets(fb60)
psi10, = wavelets(fb10)

使用尺度函数为每个滤波器组找到母小波。

sca60, = scales(fb60)
sca10, = scales(fb10)
_, idx60 = findmin(abs.(sca60 .- 1))
_, idx10 = findmin(abs.(sca10 .- 1))
m60 = psi60[idx60, :]
m10 = psi10[idx10, :]

由于 fb60 滤波器组的时间-带宽乘积更大，因此可以验证 m60 小波的包络线下的振荡比 m10 小波更多。

subplot(2, 1, 1)
plot(t, abs.(m60))
grid("on")
hold("on")
plot(t, real.(m60))
plot(t, imag.(m60))
hold("off")
xlim([-30 30])
legend("abs(m60)", "real(m60)", "imag(m60)")
title("TimeBandwidth = 60")

subplot(2, 1, 2)
plot(t, abs.(m10))
grid("on")
hold("on")
plot(t, real.(m10))
plot(t, imag.(m10))
hold("off")
xlim([-30 30])
legend("abs(m10)", "real(m10)", "imag(m10)")
title("TimeBandwidth = 10")

对齐 m60 和 m10 幅频响应的峰值。确认 m60 小波的频率响应比 m10 小波的频率响应窄。

cf60 = centerFrequencies(fb60)
cf10 = centerFrequencies(fb10)

m60cFreq = cf60[idx60]
m10cFreq = cf10[idx10]

freqShift = 2 * pi * (m60cFreq - m10cFreq);
x10 = m10 .* exp.(im .* freqShift * collect((-sigLen / 2):(sigLen / 2 - 1)));

figure()
plot(abs.(ty_fft(m60)))
hold("on")
plot(abs.(ty_fft(x10)))
grid("on")
legend("Time-bandwidth = 60", "Time-bandwidth = 10")
title("Magnitude Frequency Responses")
在多个时间序列上使用 CWT 滤波器组

本例展示了使用 CWT 滤波器组如何提高多个时间序列的 CWT 计算效率。

创建一个 100 x 1024 的矩阵 x。创建一个适合 1024 个采样点信号的 CWT 滤波器组。

using TyWavelet
using BenchmarkTools
using TyMath

rng = MT19937ar(1234)
x = randn(rng, 100, 1024)
fb = cwtfilterbank()

使用默认设置 cwt 获取 1024 个采样信号的 CWT。创建一个可包含 100 个信号的 CWT 系数的三维数组，每个信号有 1024 个采样点。

cfs, = cwt(x[1, :])
res = zeros(ComplexF64, 100, size(cfs, 1), size(cfs, 2))

使用 cwt 函数计算矩阵 x 每一行的 CWT。显示经过的时间。

@time begin
    for k in 1:100
        res[k, :, :], = cwt(x[k, :])
    end
end

1.484673 seconds (12.40 M allocations: 2.245 GiB, 14.26% gc time)

现在，使用滤波器组的 wt 对象函数求出 x 每一行的 CWT。

@time begin
    for k in 1:100
        res[k, :, :], = wt(fb, x[k, :])
    end
end


0.900339 seconds (12.34 M allocations: 1.092 GiB, 15.75% gc time)
在子图中绘制 CWT 尺度图

本示例展示了如何在数字副图中绘制 CWT 尺度图。

加载语音样本。数据采样频率为 7418 Hz。绘制默认的 CWT 尺度图。

using TyWavelet
using TyPlot
using TyBase

pkg_dir = pkgdir(TyWavelet)
source_path = pkg_dir * "/examples/Resources/mtlb.mat"
y = load(source_path)
cwt(mtlb, Fs; plotfig=true)

获取信号的连续小波变换和 CWT 频率。

cfs, frq, = cwt(mtlb, Fs)

cwt 函数在示波器中设置时间轴和频率轴。创建一个代表采样时间的向量。

tms = (0:(length(mtlb) - 1)) / Fs

在新图中，将原始信号绘制在上子图中，将尺度图绘制在下子图中。用对数标度绘制频率。

figure()
subplot(2, 1, 1)
plot(tms, mtlb)
axis("tight")
title("Signal and Scalogram")
xlabel("Time (s)")
ylabel("Amplitude")
subplot(2, 1, 2)
img = pcolor(tms, frq, abs.(cfs))
img.set_edgecolor("flat")
ax = gca()
ax.set_ylim([minimum(frq), maximum(frq)])
ax.set_xlim([tms[1], tms[end]])
ax.set_yscale("log")
xlabel("Time (s)")
ylabel("Frequency (Hz)")
ax = gca()

# 注意

  • 第一次使用滤波器组获取信号的 CWT 时，小波滤波器被构造为具有与信号相同的数据类型。将相同的滤波器组应用于不同数据类型的信号时，会生成警告消息。更改数据类型需要重新设计或更改过滤器组的精度。为了获得最佳性能，请使用一致的数据类型。
  • 当执行多个 CWT 时，例如在 for 循环中，建议的工作流程是首先创建一个 cwtfilterbank 对象，然后使用 wt 对象函数。此工作流程最大限度地减少了开销并最大限度地提高了性能。

# 参考文献

[1] Lilly, J. M., and S. C. Olhede. "Generalized Morse Wavelets as a Superfamily of Analytic Wavelets." IEEE Transactions on Signal Processing. Vol. 60, No. 11, 2012, pp. 6036–6041.

[2] Lilly, J. M., and S. C. Olhede. "Higher-Order Properties of Analytic Wavelets." IEEE Transactions on Signal Processing. Vol. 57, No. 1, 2009, pp. 146–160.

[3] Lilly, J. M. jLab: A data analysis package for MATLAB®, version 1.6.2. 2016. http://www.jmlilly.net/jmlsoft.html.

[4] Lilly, J. M. "Element analysis: a wavelet-based method for analysing time-localized events in noisy time series." Proceedings of the Royal Society A. Volume 473: 20160776, 2017, pp. 1–28. dx.doi.org/10.1098/rspa.2016.0776.

# 另请参阅

cwt | cwtfreqbounds | icwt