Matlab信号处理：连续小波变换

小波变换是信号时频分析中浓墨重彩的一笔，本文将介绍连续小波变换（Continuous Wavelet Transform，CWT），对比短时傅里叶变换(STFT)，CWT有更多的优势、更加灵活。区别于短时傅里叶变换的正弦基函数，连续小波变换采用小波基函数，通过调整小波基函数的尺度因子和时间平移因子，能分析信号在不同频率尺度和时间位置上特性。

本文的流程如下图所示，模拟了一个啁啾信号（一个线性调频信号），并借助matlab的连续小波变换函数cwt，分析了其时频特性，绘制了时频图，对比了短时傅里叶变换的结果，来表明两个方法之间的优劣。

同时，使用matlab的mesh函数、surf函数、imagesc函数、image函数、pcolor函数和imshow函数分别绘制了连续小波变换的时频图，说明不同函数之间的差异，方便大家日后作图时选择合理的绘图函数。

1.具体案例

采用短时傅里叶变换中的啁啾信号为例，啁啾(英文中为chirp)信号一般是指调频信号，即随着时间的变化，信号的主频在不断发生变化（或增加、或减少）。当这种信号当做音频输出时，听起来会像鸟的唧唧声，所以叫啁啾。 

利用matlab的chirp函数生成了一个1s内从100Hz到5000Hz的线性调频信号（啁啾信号），采样频率为24000Hz，具体如下：

从啁啾信号的局部细节图能发现，随着时间的增长，信号波形越来越密集，即信号的频率逐渐增大。在啁啾信号频谱中，低频占据信号的主频，无法发现信号频率的时变特征。

采用连续小波变换（MATLAB中cwt函数，本文采用的morse小波，其他参数均采用cwt的默认设置），获得啁啾信号的时频分析结果，如下图所示：

从上图中能发现，啁啾信号的频率从100Hz到5000Hz线性增长，这表明CWT能较好地分析该啁啾信号，对比FFT的结果，表明了CWT在信号时频分析方面的优势。

此外，在上图中CWT所表征的啁啾信号随着时间的增长，其频域的主频带宽越来越大，这似乎与啁啾信号的真实频域状态（主频的带宽不变，主频随时间线性增加）不符，这主要是因为CWT选择的小波尺度并不是最佳的。在CWT中，低频信号有更高的频率分辨率，而高频信号有更高的时间分辨率，对应的高频信号的频率分辨率就降低了，不如低频信号。简单说，随着频率的增大，其相邻频率间的间隔增大。所以，虽然啁啾信号的频谱带宽是不变的，但是在CWT的分析结果中确实逐渐增大的。如果有必要，可以自定义滤波器组或者采用STFT方法。采用短时傅里叶变换获得时频图如下：

在STFT的结果中，随着时间的增长，啁啾信号的频率逐渐增大，并且主频的带宽没发生变化，啁啾信号的时频刻画更加精细和准确，优于CWT。

CWT更适合非平稳信号的分析，它对低频信号的频率分辨率高，高频信号时间频率高，这种不同尺度的时频表征特性是我们需要的。低频信号其时间变化往往较慢，其频率多具有物理含义，如可能是转频、啮合频率、故障频率等，所以我们需要关注信号中是否有特定的低频成分，要求低频的频率分辨率高。而高频成分往往为一些杂波信号和高频噪声，这是我们需要剔除的，它们没有物理含义，所以对其频率的分辨率要求不高，多关注其出现的时刻和能量大小。（本人想法，有异议欢迎评论！）。

连续小波变换和短时傅里叶变换的绘图手段很多，结果看起来五花八门，这里简要展示一下目前在MATLAB中CWT能展示所有情况(STFT同比类似)：

mesh函数：

surf函数：

imagesc函数：

2.说明

MATLAB中适用于绘制来连续小波变换的图有影像图和三维曲面图，具体如下：

为此，将不同函数的连续小波变换结果进行了展示，其中heatmap、surfc、surfl、mehsc、meshz、waterfall、ribbon、contour3就不展示了，这些函数的具体差异，我会单独出一个表格进行简述。

上面展示了mesh函数、surf函数、imagesc函数、image函数、pcolor函数和imshow函数绘制的连续小波变换结果。其中，mesh函数、surf函数属于matlab的三维曲面图函数，imagesc函数、image函数、pcolor函数和imshow函数绘制属于matlab的影像图函数。这里需要注意：mesh函数、surf函数、pcolor函数能根据所提的X、Y坐标轴向量进行绘图，实现数值的一一对应。

而imagesc函数、image函数仅使用X、Y坐标向量的首、尾元素，按照(x(2)-x(1))/(size(C,2)-1)的方式（等间隔生成坐标轴向量，是等差数列）生成坐标向量，这导致imagesc函数、image函数在绘制连续小波变换是频率轴图像和Y轴数据不对应，请各位注意！！！为保持一致，请对y轴进行对数尺度缩放（参考mathwork的官方例子发现的，对应代码为set(gca,'yscale','log');）

mesh函数是绘制网格曲面图，即仅对三维图网格的边进行上色，概念可对比surf函数。surf函数是绘制曲面图，即对三维图网格的面进行上色，概念可对比mesh函数。下图为我放大的两个函数的绘图细节：

上图（左图为mesh函数结果，右图为surf函数结果）中，能明显发现，左图的1s处4000Hz-10000Hz的空心网格边线，而在右图的对应位置为实心曲面，这是两个函数的最大差异。

image函数是matlab的图像函数，根据二维矩阵绘制图像，它不会将矩阵最大值和最小值分别缩放到图像的颜色栏的两个极端，概念可对比imagesc函数。imagesc函数是matlab的图像函数，它会根据二维矩阵的最大值和最小值分别缩放到图像的颜色栏的两个极端，概念可对比image函数。pcolor函数绘制matlab的伪彩图，它和surf函数和mesh函数在view(0,90)视角的结果相似，但是pcolor的图不具有三维结构。值得注意的是，伪彩图上各个点Z值为0。在绘制连续小波变换结果时，记得将边颜色设置为none，否则图片可能为全黑。

imshow函数显示图像，区别于image函数、imagesc函数、pcolor函数，这三个函数会把矩阵的整体大小缩放到图窗中，这导致图横、纵轴的比例和矩阵维度的比例不一致。在本文中，连续小波变换获得矩阵大小为117*24000，但是采用image函数、imagesc函数、pcolor函数得到的图像，其X轴和Y轴似乎一样长，这与矩阵的维度时冲突的，这是因为这些函数进行了缩放，尽可能把矩阵对应的图像放到图窗中，而imshow函数不会，它得到的图像对应矩阵的维度，所以它的结果是细长的，它一般用于展示jpg、png等图片。

3.具体代码

代码主要分为两个：main.m（主函数）、frequ_am_phase.m（幅值谱和相位谱计算函数）

调用形式：[freq,P1,Theta]=frequ_am_phase(y,fs,tol)

输入：信号y，矩阵，行X列=单个信号的采样索引X信号数，比如信号的大小为8192X12，表示一个有12个信号的数据矩阵，每个信号长度为8192。注意，如果仅有一个信号，则y应该是一个列向量。同时，y的行数尽量为偶数，奇数的话会引起程序索引的警告。

输出：freq表示幅值谱的横轴，向量，HzX1；P1表示幅值谱的纵轴，矩阵，单个信号的采样索引X信号数，类比信号y；Theta表示相位谱的纵轴，矩阵，单个信号采样索引X信号数，类比信号y。

主函数main.m代码如下：

%% 信号处理——连续小波变换

clc
clear all
close all
% 定义时间向量和啁啾信号参数
fs = 24000; % 采样率
T = 1/fs:1/fs:1; % 1秒长的时间向量
f0 = 100; % 起始频率100Hz
f1 = 5000; % 结束频率5000Hz
% 生成啁啾信号
y = chirp(T, f0, 1, f1);
% 播放啁啾信号
sound(y, fs);
filename = 'chirpSignal.mp3';
% 使用audiowrite函数保存为MP3格式
audiowrite(filename, y, fs, 'Quality', 95);
%
[freq1,P1,~]=frequ_am_phase(y',fs);
figure;
subplot(311);plot(T,y,'b-');xlabel('Time/s');ylabel('Amplitude/g');axis tight;
title("啁啾信号的时域波形")
subplot(312);plot(T(1:5000),y(1:5000),'b-');xlabel('Time/s');ylabel('Amplitude/g');axis tight;
title("时域波形的局部细节")
subplot(313);plot(freq1,P1,'b-');xlabel('Frequency/Hz');ylabel('Amplitude/g');
title("啁啾信号的频域波形")
%% 连续小波变换和短时傅里叶变换结果对比
%短时傅里叶变换
[stft_result,f,t] = stft(y,fs,Window=rectwin(256),OverlapLength=128,FFTLength=256,FrequencyRange="onesided");
figure;
s1=surf(t,f,abs(stft_result),'EdgeColor','none');colorbar;
xlabel('Time/s');ylabel('Frequency/Hz');zlabel('Amplitude/g');axis tight;
title("matlab的stft函数获得时频图")
view(0,90)
%连续小波变换
[w,fre] = cwt(y,'morse',fs);
figure
s2=surf(T,fre,abs(w),'EdgeColor','none');colorbar;
xlabel('Time/s');ylabel('Frequency/Hz');zlabel('Amplitude/g');axis tight;
title("matlab的stft函数获得时频图")
view(0,90)
%% 不同方法绘图
[w,fre] = cwt(y,'morse',fs);
figure
mesh(T,fre,abs(w));colorbar;
xlabel('Time/s');ylabel('Frequency/Hz');zlabel('Amplitude/g');axis tight;
title("matlab的cwt函数获得时频图")
figure;
surf(T,fre,abs(w),'EdgeColor','none');colorbar;
xlabel('Time/s');ylabel('Frequency/Hz');zlabel('Amplitude/g');axis tight;
title("matlab的cwt函数获得时频图")
figure;
imagesc(T,fre,abs(w));axis tight;axis xy; %  连续小波图像           
set(gca,'yscale','log');            
xlabel('Time/s');ylabel('Frequency/Hz');
colorbar;
figure;
image(T,fre,abs(w));axis tight;axis xy; %  连续小波图像        
set(gca,'yscale','log');  
xlabel('Time/s');ylabel('Frequency/Hz');
colorbar;
figure;
s=pcolor(T,fre,abs(w));colorbar;
s.EdgeColor='none';
xlabel('Time/s');ylabel('Frequency/Hz');axis tight;
title("matlab的cwt函数获得时频图")
figure;
imshow(abs(w));

幅值谱和相位谱计算函数frequ_am_phase.m代码如下：

function [freq,P1,Theta]=frequ_am_phase(y,fs,tol)
%% 绘制信号频域的幅值谱和相位谱
%% 参数解释： 
%     y: 表示输入信号，它可以为一个矩阵，行X列,具体为单个信号的采样索引X信号数
%        比如y的大小为8192X12，表示一个有12个信号的数据矩阵，每个信号长度为8192
%        注意，如果仅有一个信号，则y应该是一个列向量
%        同时，y的行数尽量为偶数，奇数的话会引起程序索引的警告
%     fs:表示采样频率
%     tol:相位阈值参数
%     freq:表示幅值谱的横轴
%     P1:表示幅值谱的纵轴
%     Theta:表示相位谱的纵轴

if nargin==2
    tol=1e-6;  %计算误差的默认阈值
end

L=size(y,1);         % 信号长度
% Y=fft(y,2^nextpow2(L));          % FFT 快速傅里叶变换
Y=fft(y,L);          % FFT 快速傅里叶变换
freq=(0:L/2)*fs/L;   % 设置频率刻度  横轴Hz
%幅值谱
P2 = abs(Y/L);
P1 = P2(1:L/2+1,:);
P1(2:end-1,:) = 2*P1(2:end-1,:);  %纵轴 幅值

%相位谱
P2(2:end-1,:)=2*P2(2:end-1,:);
for i=1:size(Y,2)
    Y(P2(:,i)<tol,i) = 0;
    theta(:,i) = angle(Y(:,i))/pi;
end
Theta=theta(1:L/2+1,:);
end

综上所述，本文模拟线性调频特性的啁啾信号，并利用连续小波变换（CWT）和短时傅里叶变换(STFT)分析了它频率的时变特性；使用matlab的mesh函数、surf函数、imagesc函数、image函数、pcolor函数和imshow函数分别绘制了连续小波变换的时频图，说明不同函数之间的差异。