小波变换1-为什么要进行小波变换

文章探讨了信号处理中频率信息的重要性,通过傅里叶变换展示频率成分的检测,并指出在处理非平稳信号时,小波变换提供时频表示的优势,尤其在需要频率分量精确时间定位的场景,如脑电图分析中事件相关电位的潜伏期测量。
摘要由CSDN通过智能技术生成

提示:文章写完后,目录可以自动生成,如何生成可参考右边的帮助文档


一、变换是什么

首先,为什么我们需要变换,或者到底什么是变换?
对信号进行数学变换,以便从该信号中获取原始信号中不易获得的进一步信息。我们假设时域信号是原始信号,并将已通过任何可用的数学转换“转换”的信号假定为处理信号。
实际上,大多数信号都是原始格式的时域信号。也就是说,无论该信号测量什么,都是时间的函数。换句话说,当我们绘制信号时,其中一个轴是时间(自变量),另一个轴(因变量)通常是振幅。当我们绘制时域信号时,我们得到了信号的时域表示。对于大多数信号处理相关应用来说,这种表示并不总是信号的最佳表示。在许多情况下,最显著的信息隐藏在信号的频率成分中。信号的频谱基本上是该信号的频率分量的组成(频谱分量)。信号的频谱显示了信号中存在的频率。
直觉上,我们都知道频率与某物的速率变化有关。如果某物(数学或物理变量,在技术上是正确的术语)快速变化,我们说它是高频的,如果这个变量不会快速变化,即它平滑地变化,我们说它是低频的。如果这个变量根本没有变化,那么我们说它的频率为零,或者没有频率。例如,日报的出版频率高于月刊(出版频率更高)。
频率以周期/秒为单位,或者更常见的名称为“赫兹”。现在,请看下图,我们在matlab里产生不同频率的正弦波。第一个是 3 Hz 的正弦波,第二个是 10 Hz,第三个是 50 Hz。

%生成3Hz正弦信号
fs = 1000; %采祥率为1000HZ
t = 0:1/fs:1; %时间向量
f1 = 3; %信号频率为3Hz
x = sin(2*pi*f1*t); %生成正弦信号
%生成10Hz正弦信号
f2 = 10;%信号频率为10Hz
y = sin(2*pi*f2*t); %生成正弦信号
%生成50Hz正弦信号
f3 = 50;%信号频率为50Hz
z = sin(2*pi*f3*t); %生成正弦信号

% 创建第一个子图
subplot(3, 1, 1);
plot(t, x);
title('3hz');
xlabel('时间/s');
ylabel('幅值');
% 创建第二个子图
subplot(3, 1, 2);
plot(t, y);
title('10hz');
xlabel('时间/s');
ylabel('幅值');
% 创建第三个子图
subplot(3, 1, 3);
plot(t, z);
title('50hz');
xlabel('时间/s');
ylabel('幅值');

图1

图1

那么我们如何测量频率,或者我们如何找到信号的频率成分呢?答案是傅里叶变换(FT)。如果采用时域信号的FT,则得到该信号的频率-幅度表示。换句话说,我们现在有一个图,一个轴是频率,另一个轴是振幅。该图告诉我们信号中存在多少每个频率。
频率轴从零开始,一直到无穷大。对于每个频率,我们都有一个振幅值。例如,如果我们取第三个图,我们将有一个 50 Hz 的尖峰,而其他地方则没有,因为该信号只有 50 Hz 的频率分量。然而,没有其他信号具有如此简单的 FT。对于大多数实际用途,信号包含多个频率分量。下面显示了 50 Hz 信号的 FT:

% 定义正弦信号的参数
Fs = 1000;          % 采样频率
T = 1/Fs;            % 采样周期
t = 0:T:1;           % 时间向量
f = 50;               % 正弦信号的频率

% 生成正弦信号
x = sin(2*pi*f*t);

% 计算傅里叶变换
N = length(x);       % 信号长度
fft_result = fft(x);
frequencies = linspace(0, Fs, N);

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

图2

图2

请注意,上图中给出了两个图。底部仅绘制顶部的前半部分。实值信号的频谱始终是对称的。上图说明了这一点。但是,由于对称部分恰好是第一部分的镜像,因此它不提供其他信息,因此通常不显示此对称的第二部分。

二、为什么我们需要频率信息?

很多时候,在时域中不容易看到的信息可以在频域中看到。
让我们以生物信号为例。假设我们正在查看心电图信号(心电图,心脏电活动的图形记录)。心脏病专家熟知健康心电图信号的典型形状。任何与该形状的显着偏差通常被认为是病理状况的症状。
然而,这种病理状况在原始时域信号中可能并不总是很明显。心脏病专家通常使用记录在条形图上的时域心电图信号来分析心电图信号。最近,新的计算机化心电图记录仪/分析仪也利用频率信息来确定是否存在病理状况。当分析信号的频率成分时,有时可以更容易地诊断病理状况。
尽管FT可能是最常用的变换(尤其是在电气工程中),但它并不是唯一的变换。工程师和数学家经常使用许多其他变换。希尔伯特变换、短时傅里叶变换、当然还有我们的小波变换,只是工程师和数学家可以使用的大量变换中的一小部分。每种变换技术都有自己的应用领域,各有优缺点,小波变换(WT)也不例外。
为了更好地理解对 WT 的需求,让我们更仔细地了解 FT。FT(以及 WT)是一种可逆变换,也就是说,它允许在原始信号和处理(转换)信号之间来回切换。但是,在任何给定时间,只有其中任何一个可用。也就是说,时域信号中没有可用的频率信息,傅里叶变换信号中没有时间信息。自然而然地想到的问题是,是否有必要同时拥有时间和频率信息?
答案取决于特定的应用,以及手头信号的性质。回想一下,FT 给出了信号的频率信息,这意味着它告诉我们信号中存在多少每个频率,但它并没有告诉我们这些频率分量何时存在。当信号处于所谓的静止状态时,不需要此信息。
让我们更仔细地研究一下这个平稳性概念,因为它在信号分析中至关重要。频率成分不随时间变化的信号称为静止信号。换句话说,静止信号的频率成分不会随时间变化。在这种情况下,人们不需要知道频率分量在什么时间存在,因为所有频率分量在任何时候都存在。
例如,以下信号:
u = cos(2π10t)+cos(2π25t)+cos(2π50t)+cos(2π100t);
是一个静止信号,因为它在任何给定时间时刻的频率分别为 10、25、50 和 100 Hz。该信号如下图所示:

fs = 1000; %采祥率为1000HZ
t = 0:1/fs:0.5; %时间向量
u = cos(2*pi*10*t)+cos(2*pi*25*t)+cos(2*pi*50*t)+cos(2*pi*100*t);
plot(t, u);
title('10hz');
xlabel('时间/s');
ylabel('幅值');

在这里插入图片描述

图3

以下是它的FT:
在这里插入图片描述
图4

下图是上图的缩放版本,仅显示我们感兴趣的频率范围。请注意对应于频率 10、25、50 和 100 Hz 的四个频谱分量。
让我们看另一个例子。下图是四个不同的时间间隔绘制具有四个不同频率分量的信号,因此是 非平稳信号。间隔 0 到 300 ms 具有 100 Hz 正弦曲线,间隔 300 至 600 ms 具有 50 Hz 正弦曲线,间隔 600 至 800 ms 具有 25 Hz 正弦曲线,最后间隔 800 至 1000 ms 具有 10 Hz 正弦曲线。

clc;clear
%%
% 定义时间间隔
t1 = linspace(0, 0.3, 300);   % 0 to 300 ms
t2 = linspace(0.3, 0.6, 300); % 300 to 600 ms
t3 = linspace(0.6, 0.8, 200); % 600 to 800 ms
t4 = linspace(0.8, 1.0, 200); % 800 to 1000 ms

% 定义频率分量
frequencies = [100, 50, 25, 10];

% 生成信号
signal1 = sin(2*pi*frequencies(1)*t1);
signal2 = sin(2*pi*frequencies(2)*t2);
signal3 = sin(2*pi*frequencies(3)*t3);
signal4 = sin(2*pi*frequencies(4)*t4);

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

% 绘制信号
figure;
plot(linspace(0, 1, numel(full_signal)), full_signal);
title('非平稳信号');
xlabel('时间 (s)');
ylabel('幅度');

在这里插入图片描述

图5

以下是他的FT:
在这里插入图片描述
图6

这时不用担心小小的涟漪;它们是由于从一个频率分量到另一个频率分量的突然变化,这在本文中没有意义。请注意,高频分量的振幅高于低频分量的振幅。这是因为较高频率的分量(每个 300 毫秒)比低频分量(每个 200 毫秒)持续时间更长。(振幅的确切值并不重要)。除了那些涟漪,一切似乎都是正确的。FT有四个峰值,对应于四个具有合理振幅的频率.。
不完全错,但也不完全正确…原因如下:对于上图中绘制的第一个信号,请考虑以下问题:这些频率分量在什么时刻(或时间间隔)出现?
答:任何时候!请记住,在静止信号中,信号中存在的所有频率分量都存在于信号的整个持续时间内。始终有 10 Hz,始终有 50 Hz,始终有 100 Hz。
现在,考虑图5中 非平稳信号的相同问题。对于图5中的信号,我们知道在第一个区间中,我们的频率分量最高,在最后一个区间中,我们的频率分量最低。因此,对于这些信号,频率分量不会一直出现!比较图4和图6。这两个光谱之间的相似性应该是显而易见的。它们都显示出完全相同频率的四个频谱分量,即 10、25、50 和 100 Hz。除了纹波和振幅差异(总是可以归一化)之外,两个频谱几乎完全相同,尽管相应的时域信号甚至彼此不接近。两个信号都涉及相同的频率分量,但第一个信号始终具有这些频率,第二个信号具有不同间隔的这些频率。那么,为什么两个完全不同的信号的频谱看起来非常相似呢?回想一下,FT给出了信号的频谱内容,但它没有给出有关这些频谱分量在时间上出现的位置的信息。因此,FT不是一种适用于非平稳信号的技术,但有一个例外:
FT可以用于非平稳信号,如果我们只对信号中存在哪些频谱分量感兴趣,而不感兴趣这些分量发生在哪里。但是,如果需要这些信息,如果我们想知道,在什么时间(间隔)出现什么光谱分量,那么傅里叶变换就不是正确的变换。
出于实际目的,很难进行分离,因为有很多实用的静止信号,以及非静止信号。例如,几乎所有的生物信号都是非平稳的。一些最著名的是心电图(心脏电活动,心电图),EEG(大脑电活动,脑电图)和EMG(肌肉电活动,肌电图)。
再次请注意, FT给出了信号中存在的频率分量(频谱分量)。仅此而已。
当需要频谱分量的 时间定位时,需要进行变换,给出信号的时频表示。

三、解决方案:小波变换

小波变换是这种类型的变换。它提供了时频表示。
通常,在任何时刻出现的特定光谱分量都会引起人们的特别兴趣。在这些情况下,了解这些特定光谱分量出现的时间间隔可能非常有益。例如,在脑电图中,事件相关电位的潜伏期特别令人感兴趣(事件相关电位是大脑对特定刺激(如闪光灯)的反应,这种反应的潜伏期是刺激开始和反应之间经过的时间量)。
小波变换能够同时提供时间和频率信息,从而给出信号的时频表示。小波变换的工作原理完全是一个不同的有趣故事,应该在短时傅里叶变换 (STFT) 之后进行解释。WT 是作为 STFT 的替代品而开发的。目前只需说 WT 的开发就是为了克服 STFT 的一些与分辨率相关的问题,如第二部分所述。
我们从各种高通和低通滤波器传递时域信号,这些滤波器滤除信号的高频或低频部分。每次从信号中移除与某些频率相对应的信号的某些部分时,都会重复此过程。
假设我们有一个频率高达 1000 Hz 的信号。在第一阶段,我们将信号分成两部分,方法是将信号从高通滤波器低通滤波器(滤波器应满足某些特定条件,即所谓的允许条件)传递,从而产生同一信号的两个不同版本:对应于 0-500 Hz(低通部分)和 500-1000 Hz(高通部分)的信号部分。然后,我们取其中一部分(通常是低通部分)或两者兼而有之,并再次做同样的事情。此操作称为分解。
假设我们已经采用了低通部分,我们现在有 3 组数据,每组对应于 0-250 Hz、250-500 Hz、500-1000 Hz 频率的相同信号。无法知道信号在时频平面中某个点的频率和时间信息。换句话说:我们无法知道在任何给定时间时刻存在什么光谱成分。我们能做的最好的事情就是研究在任何给定的时间间隔内存在哪些光谱成分。这是一个分辨率问题,也是研究人员从STFT转向WT的主要原因。STFT 始终提供固定分辨率,而 WT 给出可变分辨率
较高频率在时间上分辨得更好,而较低频率在频率上分辨得更好。这意味着,与低频分量相比,某个高频分量可以更好地及时定位(相对误差更小)。相反,与高频分量相比,低频分量可以更好地定位频率。

您可以使用MATLAB中的小波变换函数来对MIT心电信号进行去噪。首先,您需要加载MIT心电信号数据,并应用小波变换。 以下是一个简单的示例代码,演示如何使用MATLAB进行小波变换去噪: ```matlab % 加载MIT心电信号数据 load mit_ecg_data.mat; % 载入小波库 load wavelets.mat; % 选择一个小波函数 waveletName = 'db4'; % 小波变换去噪参数 level = 6; % 小波分解的层数 threshold = 0.4; % 阈值 % 对每个心电信号应用小波变换去噪 denoised_signals = zeros(size(ecg_signals)); for i = 1:size(ecg_signals, 1) % 对当前心电信号应用小波变换 [C, L] = wavedec(ecg_signals(i,:), level, waveletName); % 计算软阈值 thr = threshold*sqrt(2*log(length(ecg_signals(i,:)))); % 应用软阈值 C_den = wthresh(C, 's', thr); % 重构去噪后的心电信号 denoised_signals(i,:) = waverec(C_den, L, waveletName); end % 绘制原始和去噪后的心电信号 figure; subplot(2,1,1); plot(ecg_signals(1,:)); title('原始心电信号'); subplot(2,1,2); plot(denoised_signals(1,:)); title('去噪后的心电信号'); ``` 请确保您已经将MIT心电信号数据保存在名为`mit_ecg_data.mat`的MAT文件中,并将小波库保存为名为`wavelets.mat`的MAT文件。 注意:此示例中使用了db4小波函数,您可以根据需要选择其他小波函数。另外,阈值的选择也可能需要根据具体情况进行调整。 希望这可以帮助您进行MIT心电信号的小波变换去噪。如果您有其他问题,请随时提问。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值