谨以此文，献给年少轻狂的自己

为什么需要小波变换

我们目前工程实践中所分析的信号大多是平稳信号，若忽略高斯噪声，则基本可以看作严格平稳信号的信号。但是我们难免会遇到需要进行时频分析的信号，这时候就需要采用一些时频分析方法

对于非平稳信号，传统的傅里叶变换是没有办法进行很好的分析的，需要采用特殊的方法，比如短时傅里叶变换（STFT）

短时傅里叶变换是通过添加一个固定的采样窗口来进行采样，再进行FFT，最后将结果整合为时频分析结果

下面有个例子

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
%% 预处理
clc,clf,clear;
%% 参数设置
NPT = 4000; % 信号点数
% 分段频率
f1 = 100;
f2 = 150;
f3 = 200;
f4 = 300;
sample_rate = 1000; % 采样率
t = 0:1 / sample_rate:(NPT - 1) / sample_rate;
% 数据生成
raw_signal = zeros(1,NPT);
for i = 1:1000
    raw_signal(i) = sin(2 * pi * f1 * t(i));
end
for i = 1001:2000
    raw_signal(i) = sin(2 * pi * f2 * t(i));
end
for i = 2001:3000
    raw_signal(i) = sin(2 * pi * f3 * t(i));
end
for i = 3001:4000
    raw_signal(i) = sin(2 * pi * f4 * t(i));
end
raw_signal = raw_signal + rand(1,4000);
figure('Name',"波形预览");
plot(t,raw_signal,'b');

这是一个频率和时间有关的，且混杂了高斯噪声的正弦信号的时域图

下面我们对这个信号的整个时域采样结果进行STFT分析

1
2
3
4
5
6
7
%% 进行STFT
result_stft = stft(raw_signal);
f_stft = -sample_rate / 2:sample_rate / (length(result_stft(:,1)) - 1):sample_rate / 2;
figure('Name',"STFT分析");
waterfall(0:length(result_stft(1,:)) - 1,f_stft,abs(result_stft).^2);
colorbar;
title("STFT分析"),xlabel("时间"),ylabel("频率/Hz");

从瀑布图当中可以非常明显的看出来时频关系

下面我们再对所有的采样结果进行FFT分析

1
2
3
4
5
6
%% 进行FFT(全局)
result_fft = fft(raw_signal);
f_fft = -sample_rate / 2:sample_rate / (length(result_fft) - 1):sample_rate / 2;
figure('Name',"FFT分析(全局)");
plot(f_fft,abs(fftshift(result_fft)),'b');
title("FFT分析(全局)"),xlabel("频率"),ylabel("幅值");

可以看到FFT确实也能得到频率信息，但是会丢失时间信息

但是我们知道，在实际工程当中，肯定是不能做到这种完全采集再分析的，所以更可能的情况是对任意一小段信号进行FFT

下面是随机抽取的一小段信号的FFT分析

1
2
3
4
5
6
%% 进行FFT(局部)
result_fft_ano = fft(raw_signal(1000:1350));
figure('Name',"FFT分析(局部)");
f_fft_ano = -sample_rate / 2:sample_rate / (length(result_fft_ano) - 1):sample_rate / 2;
plot(f_fft_ano,abs(fftshift(result_fft_ano)),'b');
title("FFT分析(局部)"),xlabel("频率"),ylabel("幅值");

可以看到，这次不仅丢失了时间信息，还丢失了频率信息，而这就充分的说明了时频分析方法的重要性

但是对于短时傅里叶变换而言，有一个致命的缺陷，就是窗宽不好确定

这里涉及到的就是不确定原理

即，假设ff为L2(R)L^2(R)空间的函数，它在−∞-\infty和+∞+\infty处为0，则对于任意的a∈Ra\in R和α∈R\alpha\in R，有：

Δaf⋅Δαf^≥14\Delta_af\cdot\Delta_\alpha \widehat f\ge\frac{1}{4}

意思就是ff在aa处的分辨和ff在α\alpha处的分辨不可能同时任意地小

具体证明可见《小波与傅里叶分析基础》

因为如果采样窗的长度变长了，那么时间轴的分辨率必定会下降；而采样窗的长度变短的情况下，若采样率不变，则采样点数会减少，意味着频率分辨率会下降

而且STFT的采样窗的长度是固定的，不能做到灵活可调

总的来说，缺点有二：

  1. 时间和频率分辨率都固定，不能随着频率的高低实现动态可调
  2. 选择一个合适的窗宽十分困难

而小波变换可以自由的更改时间和频率的分辨率，做到灵活可控

小波变换介绍

概念介绍

在小波分析当中，有两个函数非常重要，即尺度函数ϕ\phi和小波函数ψ\psi

这两个函数产生了一组可以用于分解和重构信号的函数族

在构造函数族的时候，ϕ\phi有时候称为父小波，ψ\psi称为母小波

需要说明，小波母函数并不是一个特定的函数，而是一种函数的集合，满足了一定条件的函数均可以作为小波母函数

小波母函数需要满足的有以下几个条件：

  1. 紧支撑性

用数学语言表示就是∃a∈>0,∀∣t∣>a,ψ(t)=0\exists a\in > 0,\forall|t| > a,\psi(t) = 0，即仅在很小的一部分定义域当中不为0，其他部分均为0

紧支撑性带来的好处是具有紧支撑性的基函数，在原信号的时间轴上平移，就相当于对于原信号就行了加窗操作

  1. 波动性

即∫−∞+∞ψ(t)dt=0\int^{+\infty}_{-\infty}\psi(t)dt = 0，也就是说，在所有定义域内积分为0

  1. 容许条件

cψ=2π∫−∞∞∣ψ(f~)∣2∣f∣df=0c_\psi = 2\pi\int^{\infty}_{-\infty}\frac{|\psi(\widetilde f)|^2}{|f|}df = 0

其中，ψ(f~)\psi(\widetilde f)是小波函数傅里叶变换的共轭，这个条件使得变换可逆

  1. 正交性

这个条件也是为了变换可逆

过程分析

小波母函数的变换公式如下

ψ∗(τ,s)=1sψ(ts−τ)\psi^\ast(\tau,s) = \frac{1}{\sqrt s}\psi(\frac{t}{s} - \tau)

这体现出了小波变换的两种变换方式：平移和缩放

平移是通过τ\tau控制的，缩放是通过ss控制的

再平移之后，就需要和信号相乘了，这可以筛选出在小波母函数非零部分，频率和小波母函数相近的成分

将上述过程写成公式就是

CWTxψ(τ,s)=1s∫−∞∞x(t)ψ(ts−τ)dtCWT^\psi_x(\tau,s) = \frac{1}{\sqrt s}\int^{\infty}_{-\infty}x(t)\psi(\frac t s - \tau)dt

事实上，τ\tau和ss可以为任意数，即有无限种组合

仿真对比

1
2
3
4
5
6
7
%% CWT
cwt_result = cwt(raw_signal);
f_cwt = -sample_rate / 2:sample_rate / (length(cwt_result(:,1)) - 1):sample_rate / 2;
figure('Name',"CWT分析");
waterfall(0:length(cwt_result(1,:)) - 1,f_cwt,abs(cwt_result).^2);
colorbar;
title("CWT分析");

可以看到CWT的分辨率是优于STFT的

常用小波基介绍

我们常用的小波基函数就以下几种

小波基函数名称 优点 缺点
Haar 简单 性能一般
Daubechies(dbN) 常用于滤波器 性能和迭代次数有关
Mexican Hat(mexh) 在时间域与频率域都有很好的局部化 未知
Morlet 未知 未知
Meyer 未知 未知

下面会解释简单的Haar小波基

Haar

概念

Haar的尺度函数为

ϕ(x)={1,for 0≤x<10,otherwise\phi(x) = \begin{cases}1,&\text{for $0\le x< 1$}\\\\0,&\text{otherwise}\end{cases}

小波函数为

ψ(x)={1,for 0≤x<12−1,for 12≤x<10,otherwise\psi(x) = \begin{cases}1,&\text{for $0\le x< \frac 12$}\\\\-1,&\text{for $\frac 12\le x<1$}\\\\0,&\text{otherwise}\end{cases}

令V0V_0是所有形如

∑k∈Zakϕ(x−k),ak∈R\sum_{k\in Z}a_k\phi(x - k),a_k\in R

的函数组成的空间，kk可以在任一有限的整数范围内取值

因为ϕ(x−k)\phi(x - k)在x=kx = k和x=k+1x = k + 1处不连续，换句话说，V0V_0是所有不连续点仅在整数集中的分段常量函数所组成的空间

当然，V0V_0的函数也可能在整数点处连续，如a1=a2a_1 = a_2的情况，和式在x=2x = 2处连续

设jj是非负整数，jj级阶梯空间表示为VjV_j，它是由函数集

{⋯ ,ϕ(2jx+1),ϕ(2jx),ϕ(2jx−1),ϕ(2jx−2),⋯ }\{\cdots,\phi(2^jx + 1),\phi(2^jx),\phi(2^jx - 1),\phi(2^jx - 2),\cdots\}

在实数域上张成的，VjV_j是紧支撑的分段常量函数空间，其间断点在下列集合中：

{⋯ ,−12j,0,12j,22j,32j,⋯ }\{\cdots,-\frac{1}{2^j},0,\frac{1}{2^j},\frac{2}{2^j},\frac{3}{2^j},\cdots\}

V0V_0中的函数是在整数集上有间断点的分段常量函数，V0V_0中任何一个函数亦属于V1V_1，依此类推，有：

V0⊂V1⊂V2⊂V3⊂⋯Vj−1⊂Vj⊂Vj+1⊂⋯V_0\subset V_1\subset V_2\subset V_3\subset\cdots V_{j - 1}\subset V_j\subset V_{j + 1}\subset\cdots

这种包含关系是严格的，VjV_j包含所有在分辨率为2−j2^{-j}下的相关信息，随着jj的增加，分辨的更精细

而Vj⊂Vj+1V_j\subset V_{j + 1}意味着随分辨率的提高，不会损失任何信息，该包含关系也说明了为什么VjV_j是以ϕ(2jx)\phi(2^jx)形式而不是以ϕ(ax)\phi(ax)形式定义的

下面不经过证明给出一个定理

函数集{2j2ϕ(2jx−k);k∈Z}是Vj的一个标准正交基函数集\{2^{\frac{j}{2}}\phi(2^jx - k);k\in Z\}是V_j的一个标准正交基

但是我们有了VjV_j的正交基之后，也只是完成了一半的任务

为了解决噪声滤除问题，我们需要一个孤立属于VjV_j但是不属于Vj−1V_{j - 1}的方法，于是引出了小波函数ψ\psi

该方法就是将VjV_j分解成Vj−1V_{j - 1}及其正交补

令WjW_j是由下列函数构成的空间

∑k∈Zakψ(2jx−k)ak∈R\sum_{k\in Z}a_k\psi(2^jx- k)\quad a_k\in R

这里假定只有有限个aka_k非零，WjW_j是VjV_j的正交补，即

Vj+1=Vj⊕WjV_{j + 1} = V_j\oplus W_j

也就是说VjV_j可以不断分解下去，直到0，最后得到以下结果

VjV_j中的任一ff可唯一地分解为以下的和式

f=wj−1+wj−2+⋯+w0+f0f = w_{j - 1} + w_{j - 2} + \cdots + w_0 + f_0

其中wl∈Wl,0≤l≤j−1w_l\in W_l,0\le l\le j - 1且f0∈V0f_0\in V_0

直观的说，wlw_l表示宽为12l+1\frac{1}{2^{l + 1}}的”尖峰“，且不能由其他宽度的尖峰的线性组合表示

我们再次不加证明地给出以下定理

L2(R)能被分解为无限个正交直和:L2(R)=V0⊕W0⊕W1⊕⋯特别地，对于每个f∈L2(R)可以唯一的写成:f=f0+∑j=0∞wj此处f0∈V0,wj∈Wj无限直和可以视为有限直和的极限，即f=f0+lim⁡N→∞wj这里的极限是指L2意义上的极限L^2(R)能被分解为无限个正交直和:\\L^2(R) = V_0\oplus W_0\oplus W_1\oplus\cdots\\特别地，对于每个f\in L^2(R)可以唯一的写成:\\f = f_0 + \sum^\infty_{j = 0}w_j\\此处f_0\in V_0,w_j\in W_j\\无限直和可以视为有限直和的极限，即\\f = f_0 + \lim_{N\rightarrow\infty}w_j\\这里的极限是指L^2意义上的极限

分解

至此，从理论上讲，噪声的滤除问题已经解决了

fjf_j可以分解为：

fj=f0+w1+w2+⋯+wj−1,w1∈W1f_j = f_0 + w_1 + w_2 + \cdots + w_{j - 1},\quad w_1\in W_1

其中wlw_l表示宽为12l+1\displaystyle\frac{1}{2^{l + 1}}的尖峰，当ll足够大时，这样的尖峰就能够表示噪声了

下面我们来实现分解算法

先用以下的阶梯函数近似原信号ff：

fj(x)=∑i∈Zaiϕ(2jx−l)f_j(x) = \sum_{i\in Z}a_i\phi(2^jx - l)

该过程就是在对信号采样，采样点为x=⋯ ,−12j,0,12j,⋯x = \cdots,-\displaystyle\frac{1}{2^j},0,\displaystyle\frac{1}{2^j},\cdots从而得到al=f(l2l),l∈Za_l = f(\frac{l}{2^l}),l\in Z

下面不加证明地给出ϕ\phi和ψ\psi之间的关系式

ϕ(2jx)=(ψ(2j−1x)+ϕ(2j−1x))/2ϕ(2jx−1)=(2j−1ϕ(x)−2j−1ψ(x))/2\phi(2^jx) = (\psi(2^{j - 1}x) + \phi(2^{j - 1}x)) / 2\\\phi(2^jx - 1) = (2^{j - 1}\phi(x) - 2^{j - 1}\psi(x)) / 2

可以由这个关系式将原信号进行分解，最后得到以下分解公式

设fj(x)=∑k∈Zakjϕ(2jx−k)∈Vj那么fj可以分解为:fj=wj−1+fj−1这里,wj−1=∑k∈Zbkj−1ψ(2j−1x−k)∈Wj−1fj−1=∑k∈Zakj−1ψ(2j−1x−k)∈Vj−1其中,bkj−1=a2kj−a2k+1j2akj−1=a2kj+a2k+1j2该分解过程可以用j−1取代j,继续分解下去,最后得到fj=wj−1+wj−2+⋯+w0+f0设\\f_j(x) = \sum_{k\in Z}a^j_k\phi(2^jx - k) \in V_j\\那么f_j可以分解为:\\f_j = w_{j - 1} +f_{j - 1}\\这里,\\w_{j - 1} = \sum_{k\in Z}b_k^{j - 1}\psi(2^{j - 1}x - k)\in W_{j - 1}\\f_{j - 1} = \sum_{k\in Z}a_k^{j - 1}\psi(2^{j - 1}x - k)\in V_{j - 1}\\其中,\\b_k^{j - 1} = \frac{a^j_{2k} - a^j_{2k + 1}}{2}\quad a^{j - 1}_k = \frac{a^{j}_{2k} + a^j_{2k + 1}}{2}\\该分解过程可以用j - 1取代j,继续分解下去,最后得到\\f_j = w_{j - 1} + w_{j - 2} + \cdots + w_0 + f_0

重构

至此，我们已经完成信号的分解

而我们最终的目的是信号滤波，所以下一步应该将信号重构回去

但是在此之前，我们需要丢弃一部分被视为噪声的WjW_j

如果是为了压缩数据，则较小幅值的WjW_j可以被丢弃，仅传输显著的分量，从而获得较大的压缩比

重构算法的证明可以见《小波与傅里叶分析基础》

下面我们不加证明地给出Haar重构公式

设f=f0+w0+w1+⋯+wj−1这里,f0(x)=∑k∈Zak0ϕ(x−k)∈V0wj′(x)=∑k∈Zbkj′ψ(2j′x−k∈Wj′)其中0≤j′<j,那么f(x)−∑l∈Zaijϕ(2jx−l)∈Vj这里的alj′是根据如下算法,由j′=1,j′=2,⋯循环确定的,alj′−{akj′−1+bkj′−1,若l=2k是偶数akj′−1−bkj′−1,若l=2k+1是奇数设\\f = f_0 + w_0 + w_1 + \cdots + w_{j - 1}\\这里,\\f_0(x) = \sum_{k\in Z}a^0_k\phi(x - k)\in V_0\quad w_{j'}(x) = \sum_{k\in Z}b^{j'}_k\psi(2^{j'}x - k\in W_{j'})\\其中0\le j'<j,那么\\f(x) - \sum_{l\in Z}a^j_i\phi(2^jx - l)\in V_j\\这里的a^{j'}_l是根据如下算法,由j' = 1,j' = 2,\cdots循环确定的,\\a^{j'}_l - \begin{cases}a^{j' - 1}_k + b^{j' - 1}_k,&\text{若$l = 2k$是偶数}\\\\a^{j' - 1}_k - b^{j' - 1}_k,&\text{若$l = 2k +1$是奇数}\end{cases}

多分辨率分析

暂略

道阻且长，To be continued⋯\cdots