心电信号傅里叶变换的理解

upbobo

已于 2023-04-10 11:14:34 修改

阅读量995

点赞数 2

分类专栏：心电信号处理文章标签： c++ python matlab

于 2023-03-30 20:24:13 首次发布

本文链接：https://blog.csdn.net/m0_58433013/article/details/129865100

版权

心电信号处理专栏收录该内容

1 篇文章 0 订阅

订阅专栏

从MIT-BIH数据库下载心电数据200m.mat和233m.mat

使用plotATM函数（需要自己编写后封装）对心电数据进行处理，处理后的每组心电数据长度N=2048，分别命名为ECGval200.mat和ECGval233.mat。下面贴上plotATM函数代码及调用plotATM函数的程序（是MATLAB程序哦）

plotATM函数代码：

function val = plotATM(Name)
infoName = strcat(Name, '.info'); %水平串联字符串
matName = strcat(Name, '.mat');
Octave = exist('OCTAVE_VERSION'); %检查变量、脚本、函数、文件夹或类的存在情况
load(matName);
fid = fopen(infoName, 'rt');
fgetl(fid); %读取文件中的行，并删除换行符
fgetl(fid);
fgetl(fid);
[freqint] = sscanf(fgetl(fid), 'Sampling frequency: %f Hz  Sampling interval: %f sec');
interval = freqint(2);
fgetl(fid);

if(Octave)
    for i = 1:size(val, 1)
       R = strsplit(fgetl(fid), char(9));
       signal{i} = R{2};
       gain(i) = str2num(R{3});
       base(i) = str2num(R{4});
       units{i} = R{5};
    end
else
    for i = 1:size(val, 1)
      [row(i), signal(i), gain(i), base(i), units(i)]=strread(fgetl(fid),'%d%s%f%f%s','delimiter','\t');
    end
end

fclose(fid);
val(val==-32768) = NaN;

for i = 1:size(val, 1)
    val(i, :) = (val(i, :) - base(i)) / gain(i);
end

x = (1:size(val, 2)) * interval;
plot(x', val');

for i = 1:length(signal)
    labels{i} = strcat(signal{i}, ' (', units{i}, ')'); 
end

legend(labels); %在坐标区上添加图例
xlabel('Time (sec)');
end

调用plotATM函数的程序：

%% 
val = plotATM('200m');
ECGval=val(1:2048);
save ECGval200.mat ECGval

%% 
val = plotATM('233m');
ECGval=val(1:2048);
save ECGval233.mat ECGval

随后绘制ECGval200.mat和ECGval233.mat的原始心电波形图，代码如下：

clear
Fs=360;N=2048;
load('ECGval200.mat');
x1=ECGval;
load('ECGval233.mat');
x2=ECGval;
figure
subplot(2,1,1);plot(x1);grid;title('原始心电信号x1.200.mat');
subplot(2,1,2);plot(x2);title('原始心电信号x2.233.mat');grid

结果图：

！切入本篇文章的正题：对ECGval200.mat和ECGval233.mat两个心电信号进行傅立叶变换

傅立叶变换可以将信号从时域转换到频域，从时域的y-t图形转换到频域的幅值-频率图形，得到信号的各个频率分量。代码：

mf1=fft(x1,N);%进行频谱变换（傅里叶变换）
mf2=fft(x2,N);
mag1=abs(mf1);%取绝对值
mag2=abs(mf2);
f1=(0:length(mf1)-1)*Fs/length(mf1);  %进行频率变换   这是频谱图的横坐标
f2=(0:length(mf2)-1)*Fs/length(mf2);

变换前x1、x2数据长度为N，傅立叶变换后mf1、mf2数据长度也为N，对mf1和mf2取绝对值，就得到了2048（前面设置的N）个正数，这就是傅立叶变换后，得到的信号在2048个频率分量下的幅度特性。

此时如果使用plot绘制mf1或mf2的图形，可以得到横坐标为0~2047、纵坐标为幅度的折线图。下图所示。

若想得到频谱图，则需将横坐标转换为对应的频率值，这就需要对横坐标0~2047进行频率变换，频率变换的过程：横坐标*采样率/采样点数。

采样率/采样点数=频谱图中相邻两点之间的频率之差。

如何理解呢？

采样率Fs，采样点数N，则完成N点采样所需时间 $T=\frac{1}{Fs}*N$ ，

因此这N个点的频率分辨率 $F=\frac{1}{T}=\frac{Fs}{N}$ ，即相邻两点之间对应的频率值相差Fs/N。

频谱图中，第1个点的横坐标是0，代表着傅里叶变换后得到的直流分量；第i个点的横坐标是 $\left ( i-1 \right )*\frac{Fs}{N}$ ，代表着傅里叶变换后得到的第i个频率分量。绘制幅度谱：

可以看到，横坐标是频率值，我们可以轻易地知道傅立叶变换后的频率分量及其幅度

傅立叶变换后得到的频谱图是关于 $\frac{Fs}{2}$ 对称的。

绘制幅度谱和相位谱的代码放下面：

figure
subplot(2,2,1);
plot(f1,mag1);axis([0,370,0,190]);grid;      %画出频谱图  
xlabel('频率(HZ)');ylabel('幅值');title('心电信号x1幅度谱'); 
subplot(2,2,2);
plot(f2,mag2);axis([0,370,0,190]);grid;      %画出频谱图  
xlabel('频率(HZ)');ylabel('幅值');title('心电信号x2幅度图'); 
subplot(2,2,3);
plot(f1,angle(mf1)/pi);axis([0 360 -1.3 1.3]);grid;
xlabel('频率(HZ)');ylabel('相位');title('心电信号x1相位谱'); 
subplot(2,2,4);
plot(f1,angle(mf2)/pi);axis([0 360 -1.3 1.3]);grid;
xlabel('频率(HZ)');ylabel('相位');title('心电信号x2相位谱');

后续会分享更多有关心电信号处理的内容，有时间了会码上来的~