小波变换2-傅里叶变换和短时傅里叶变换

一、基础知识

让我们对第一部分进行简短回顾。我们基本上需要小波变换（WT）来分析非平稳信号，即频率响应随时间变化的信号。我已经写过，傅里叶变换（FT）不适合非平稳信号，并且我已经展示了它的示例以使其更加清楚。为了快速回忆一下，我举一个例子。假设我们有两个不同的信号。还假设它们都具有相同的光谱成分，但有一个主要区别。假设其中一个信号始终具有四个频率分量，而另一个信号在不同时间具有相同的四个频率分量。两个信号的 FT 相同，如第 1 部分的示例所示。尽管这两个信号完全不同，但它们的 FT（幅度）是相同的！这显然告诉我们不能将 FT 用于非平稳信号。
但为什么会出现这种情况呢？换句话说，为什么两个信号具有相同的 FT？傅里叶变换到底是如何工作的？

二、傅里叶变换

在 19 世纪，法国数学家 J. Fourier 表明，任何周期函数都可以表示为表示为周期性复指数函数的无限和。在他发现（周期）函数的这一显着性质多年后，他的想法首先被推广到非周期函数，然后推广到周期或非周期离散时间信号。正是在这种概括之后，它成为一种非常适合计算机计算的工具。1965年，一种称为快速傅里叶变换（FFT）的新算法被开发出来，FT变得更加流行。
傅里叶变换告诉我们某个频率分量是否存在。该信息与该组件出现的时间无关。因此，在使用 FT 处理信号之前，了解信号是否平稳非常重要。
第一部分中给出的例子现在应该很清楚了。我在这里再给一下：
看下图，它显示了信号：

u = cos(2π10t)+cos(2π25t)+cos(2π50t)+cos(2π100t);

也就是说，它具有 5、10、20 和 50 Hz 四种频率分量，并且始终出现。

fs = 1000; %采祥率为1000HZ
t = 0:1/fs:0.5; %时间向量
u = cos(2*pi*5*t)+cos(2*pi*10*t)+cos(2*pi*20*t)+cos(2*pi*50*t);
plot(t, u);
title('10hz');
xlabel('时间/s');
ylabel('幅值');
% 计算傅里叶变换
N = length(u);       % 信号长度
fft_result = fft(u);
frequencies = linspace(0, Fs, N);

% 绘制傅里叶变换的幅度谱
figure;
plot(frequencies(1:120), abs(fft_result(1:120)));
title('傅里叶变换幅度谱');
xlabel('频率 (Hz)');
ylabel('幅度');

这里频率轴被切掉了，但理论上它延伸到无穷大（对于连续傅里叶变换（CFT）。实际上，这里我们计算离散傅里叶变换（DFT），在这种情况下频率轴上升（至少）两倍信号的采样频率，变换后的信号是对称的。不过，此时这并不那么重要。）

注意上图中的四个峰值，它们对应于四个不同的频率。
通过第一部分，现在我们知道我们不能（好吧，我们可以，但我们不应该）将 FT 用于非平稳信号，我们该怎么办？

三、短时傅里叶变换

FT出了什么问题？它不适用于非平稳信号。让我们这样想：我们可以假设非平稳信号的某些部分是平稳的吗？答案是肯定的。
如果可以假设信号平稳的区域太小，那么我们就从狭窄的窗口中观察该信号，该窗口足够窄，以至于从这些窗口看到的信号部分确实是平稳的，即所谓的：短时傅里叶变换（STFT）。
STFT和FT之间只有很小的区别。在 STFT 中，信号被分成足够小的段，其中信号的这些段（部分）可以假设是静止的。为此，选择窗函数“w”。该窗口的宽度必须等于信号平稳性有效的部分。
该窗口函数首先位于信号的最开头。即，窗函数位于t=0处。假设窗口的宽度是“T”。在这个时刻（t=0），窗口函数将与前 T/2 秒重叠（我假设所有时间单位都是秒）。然后将窗函数和信号相乘。通过这样做，仅选择信号的前 T/2 秒，并对窗口进行适当的加权（如果窗口是矩形，幅度为“1”，则乘积将等于信号）。然后假设该乘积只是另一个信号，需要对其进行 FT。换句话说，取该乘积的 FT，就像取任何信号的 FT 一样。
该窗口函数首先位于信号的最开头。即，窗函数位于t=0处。假设窗口的宽度是“T”。在这个时刻（t=0），窗口函数将与前 T/2 秒重叠（我假设所有时间单位都是秒）。然后将窗函数和信号相乘。通过这样做，仅选择信号的前 T/2 秒，并对窗口进行适当的加权（如果窗口是矩形，幅度为“1”，则乘积将等于信号）。然后假设该乘积只是另一个信号，需要对其进行 FT。换句话说，取该乘积的 FT，就像取任何信号的 FT 一样。
此转换的结果是信号的前 T/2 秒的 FT。如果信号的这一部分是静止的，正如所假设的那样，那么就不会有问题，并且获得的结果将是信号的前 T/2 秒的真实频率表示。
下一步，将该窗口（大约 t1 秒）移动到新位置，与信号相乘，并取乘积的 FT。遵循此过程，直到通过以“t1”秒间隔移动窗口到达信号末尾。
下图可能会帮助您更好地理解这一点：

颜色中的类高斯函数是窗函数。红色显示位于 t = t1’ 的窗口，蓝色显示 t = t2’ 的窗口，绿色显示位于 t = t3’ 的窗口。这些将对应于三个不同时间的三个不同的 FT。因此，我们将获得信号的 真实时频表示（TFR） 。
首先，由于我们的变换是时间和频率的函数（与 FT 不同，FT 仅是频率的函数），因此变换将是二维的（如果也计算幅度，则为三维）。让我们采用一个非平稳信号，下图是四个不同的时间间隔绘制具有四个不同频率分量的信号，因此是 非平稳信号。在这个信号中，有四个不同时间的频率分量。0 到 250 ms 的间隔是 300 Hz 的简单正弦波，其他 250 ms 间隔分别是 200 Hz、100 Hz 和 50 Hz 的正弦波。显然，这是一个非平稳信号。现在，让我们看看它的 STFT：。

%% 生成示例信号
Fs = 1000;  % 采样频率
t = 0:1/Fs:1;  % 时间向量，总共1秒

% 0 到 250 ms 间隔是 300 Hz 的简单正弦波
signal1 = sin(2*pi*300*t(1:round(0.25*Fs)));

% 其他 250 ms 间隔分别是 200 Hz、100 Hz 和 50 Hz 的正弦波
signal2 = sin(2*pi*200*t(round(0.25*Fs)+1:round(0.5*Fs)));
signal3 = sin(2*pi*100*t(round(0.5*Fs)+1:round(0.75*Fs)));
signal4 = sin(2*pi*50*t(round(0.75*Fs)+1:end));

% 合并信号
signal = [signal1, signal2, signal3, signal4];

% 绘制信号曲线
figure;
plot(t, signal);
title('信号的时域表示');
xlabel('时间 (s)');
ylabel('幅度');

%% 计算 STFT
window_length = round(Fs*0.05);  % 窗口长度，选择0.05秒
overlap_length = round(window_length/2);  % 重叠长度，选择一半的窗口长度
[S, F, T] = spectrogram(signal, window_length, overlap_length, [], Fs);

% 绘制STFT的三维图
figure;
surf(T, F, (abs(S)), 'EdgeColor', 'none');
title('信号的短时傅里叶变换（STFT）');
xlabel('时间 (s)');
ylabel('频率 (Hz)');
zlabel('振幅 (dB)');
colorbar;

显然，这是一个非平稳信号。现在，让我们看看它的 STFT：

有四个峰值对应于四个不同的频率分量。与 FT 不同，这四个峰值沿时间轴位于不同的时间间隔。请记住，原始信号具有位于不同时间的四个频谱分量。现在我们有了信号的真实时频表示。
我们不仅知道信号中存在哪些频率分量，而且还知道它们在时间上的位置。既然 STFT 给出了信号的 TFR，为什么我们还需要小波变换。STFT 隐含的问题在上面的例子中并不明显。
STFT 的问题在于其根源可以追溯到所谓的海森堡不确定性原理。该原理最初应用于运动粒子的动量和位置，现在也可以应用于信号的时频信息。简而言之，该原理指出，人们无法知道信号的准确时频表示，即无法知道在什么时间存在什么频谱分量。我们能知道的是某个频段存在的时间间隔，这是一个分辨率问题。 STFT 的问题与所使用的窗函数的宽度有关。
在 FT 中，频域中不存在分辨率问题，即我们确切地知道存在哪些频率；类似地，时域中不存在时间分辨率问题，因为我们知道每个时刻的信号值。FT 中提供完美频率分辨率的事实是，FT 中使用的窗口是其内核，该函数始终从负无穷大到正无穷大。现在，在 STFT 中，我们的窗口长度有限，因此它仅覆盖信号的一部分，这导致频率分辨率变差。变得更差是指，我们不再知道信号中存在的确切频率分量，但我们只知道存在的一个频段。在 FT 中，核函数使我们能够获得完美的频率分辨率，因为内核本身是一个无限长度的窗口。在 STFT 中，窗口的长度是有限的，我们不再具有完美的频率分辨率。为什么我们不像 FT 那样将 STFT 中的窗口长度设置为无限大，以获得完美的频率分辨率？好吧，如果你一直丢失信息，你基本上会得到 FT 而不是 STFT。
如果我们使用无限长度的窗口，我们得到的 FT 会给出完美的频率分辨率，但没有时间信息。此外，为了获得平稳性，我们必须有一个足够短的窗口，其中信号是平稳的。窗口越窄，时间分辨率越好，平稳性假设越好，但频率分辨率越差：
窄窗口===>时间分辨率好，频率分辨率差。
宽窗口===>频率分辨率好，时间分辨率差。
下图显示了不同宽度的四个窗口函数：

%% 计算 STFT
% 不同宽度的窗口函数
window_lengths = [round(Fs*0.025),round(Fs*0.05), round(Fs*0.1), round(Fs*0.2)];  % 不同窗口宽度

window_length = (round(Fs*0.1));  
for i = 1:length(window_lengths)
    window_length = window_lengths(i);
    overlap_length = round(window_length/2);  % 重叠长度，选择一半的窗口长度
    
    % 计算并绘制STFT

    [S, F, T]=spectrogram(signal, window_length, overlap_length, [], Fs);
    
    % 绘制STFT的三维图
    subplot(2, 2, i);
%     figure;
    surf(T, F, (abs(S)), 'EdgeColor', 'none');
    title(['STFT with Hann Window - Window Length: ' num2str(window_length) ' samples']);
    xlabel('时间 (s)');
    ylabel('频率 (Hz)');
    zlabel('振幅 (dB)');
    view(-45,45);  
    colorbar;
end

首先看看第一个最窄的窗口。我们期望 STFT 具有非常好的时间分辨率，但频率分辨率相对较差：

上图展示了这个STFT。该图是从顶部鸟瞰图显示的，并带有一个角度以便更好地解释。请注意，四个峰值在时间上彼此分开。另请注意，在频域中，每个峰值都覆盖一个频率范围，而不是单个频率值。现在让我们将窗口加宽，并查看第三个窗口（第二个窗口已在第一个示例中显示）。

与之前的情况不同，峰值在时间上彼此分离得并不好，但是，在频域中，分辨率要好得多。现在让我们进一步增加窗口的宽度，看看会发生什么：

时间分辨率会变得更差。
这些例子应该说明了 STFT 隐含的分辨率问题。任何想要使用 STFT 的人都会面临这个分辨率问题。使用什么样的窗户？窄窗口提供良好的时间分辨率，但频率分辨率较差。宽窗口提供良好的频率分辨率，但时间分辨率较差；此外，宽窗可能会违反平稳性条件。当然，问题是一劳永逸地选择窗口函数并在整个分析中使用该窗口的结果。当然，答案取决于应用：如果原始信号中的频率分量彼此很好地分离，那么我们可能会牺牲一些频率分辨率并追求良好的时间分辨率，因为频谱分量已经彼此很好地分离。然而，如果情况并非如此。
小波变换（WT）在一定程度上解决了分辨率的困境，我们将在下一部分中看到。