一、基础知识
让我们对第一部分进行简短回顾。我们基本上需要小波变换(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)在一定程度上解决了分辨率的困境,我们将在下一部分中看到。