摘要

本研究通过采用傅里叶变换方法，利用Python编程环境对帕金森病患者的步态数据进行详细分析，以探讨疾病对患者步态的影响。研究目的在于通过收集和分析步态参数（如步长、步速及步态不对称性）开发一套辅助诊断帕金森病、监测病情进展，并评估治疗效果的方法。实验在Google Colab平台进行，涵盖信号生成、手动与自动傅里叶变换及步态信息的时域和频域分析。研究结果表明，帕金森病患者与健康对照组在步态特征上存在显著差异，特别是在频域分析中明显体现。结论指出，傅里叶变换不仅能有效分析和解释复杂生物信号，其在医学诊断和生物力学研究中的应用具有重要价值。

一、引言

1.背景知识

在现代科学和工程的许多领域中，对信号进行分析与处理是不可或缺的一环。信号可以是任何形式的信息传递，如音频、视频、地震波或医学数据。为了深入理解这些信号所含的信息，科学家和工程师需要一种工具来分析其频率成分，即信号中各个频率的强度和相位。这种需求促进了傅里叶变换（Fourier Transform, FT）的发展，它是由18世纪末的法国数学家让-巴普蒂斯·傅里叶提出的。
傅里叶变换是一种数学变换，用于将信号从时间域（或空间域）转换到频率域。这种变换不仅揭示了信号的周期性和非周期性成分，而且还提供了一种分析和处理复杂信号的方法。例如，在连续傅里叶变换中，任何实际的物理信号都可以表达为无限多的正弦波和余弦波的积分，这种表示形式极大地便利了频率分析。
傅里叶变换的应用广泛而深远，涵盖了从电子工程到量子物理的多个领域。在电信领域，傅里叶变换帮助工程师分析和设计滤波器和信号编解码器。在图像处理领域，通过频率域分析，可以有效地进行图像压缩和图像增强。此外，傅里叶变换在医学成像技术中也扮演了关键角色，如磁共振成像（MRI）和计算机断层扫描（CT）。
由于傅里叶变换在处理非周期性和时间变化信号方面的独特优势，它在现代科技发展中占据了核心地位。深入理解傅里叶变换的理论和应用，不仅对学术研究具有重要意义，也对工业实践和技术创新具有指导作用。

2.实验目的和意义

本次实验的主要目的是利用Python编程语言来分析帕金森病患者的步态数据，以识别和量化帕金森病对患者步态的影响。通过收集和分析步态参数，如步长、步速、步态不对称性等，本实验旨在开发一种有效的方法来辅助诊断帕金森病，监测病情进展，并评估治疗效果。学术意义方面，首先，步态分析提供了一种非侵入性的方法来研究帕金森病患者的运动障碍，这是理解帕金森病中运动功能衰退的重要方面。此实验有助于揭示步态变化与帕金森病病理之间的关联，从而促进对该疾病更深层次的认识。通过基于傅里叶变换的步态信息分析通过将步态数据从时域转换到频域， 提取关键的频率成分和步态特征， 为步态评估和疾病诊断提供了一种有效的工具。 这种方法能够揭示步态的内在规律和异常模式， 对于理解人类行走机制和促进相关领域的研究具有重要意义。

二、实验方法

1.实验环境

Google Colab

2.实验步骤

2.1 生成信号，进行手动傅里叶变换以及内置 FFT 函数傅里叶变换

生成信号： 通过自回归过程生成信号$x\left(t\right)={\sum }_{i=1}^k\left({\alpha }_{i}x\left(t-1\right)\right)+\sin 0.05\pi t+{ϵ}_t$

# temporal weighting vector
w = np.array([-.6,.9]) # -2,-1
k = len(w)
N = 200
x = np.zeros(N)
for i in range(k,N):x[i] = sum( w*x[i-k:i] ) + np.random.randn()
# add sine wave
x += np.sin(np.linspace(0,10*np.pi,N))# plot the signal
plt.plot(x)
plt.xlabel('X 轴')
plt.ylabel('函数值')
plt.title('通过自回归过程生成信号')
plt.show()

这段代码首先定义了一个自回归(AR)模型的权重向量w，包含两个元素 -0.6 和 0.9，用来指定自回归过程中前两个时间点的权重。接着，它初始化了一个长度为 200 的信号数组 x，数组的前两个值默认为零。随后，对每个时间点从 k 开始到 N，利用前两个时间点的加权和以及一个随机噪声（正态分布）来计算当前时间点的信号值。之后，代码在这个自回归过程生成的信号上叠加了一个正弦波，其频率随着时间线性增加。最后，使用 Matplotlib 绘制了这个信号，并设置了图表的各种标签和标题

2.2 进行手动傅里叶变换以及内置 FFT 函数傅里叶变换

其中手动傅里叶变换通过将信号$x(t)$与不同频率的$e^{-j2\pi ft}$进行内积计算得到。

# time vector
t = np.arange(0,N)/N# initialize Fourier coefficients
fc = np.zeros(N,dtype=complex)for f in range(N):# create complex sine wavecsw = np.exp(-1j*2*np.pi*f*t)# dot product with signalfc[f] = np.dot(csw,x)# check against the fft
fc2 = np.fft.fft(x)# normalized frequencies vector (avoid plotting the negative frequencies)
hz = np.linspace(0,1,int(N/2+1))# and plot
plt.plot(hz,np.abs(fc[:len(hz)]),label='手动 FT')
plt.plot(hz,np.abs(fc2[:len(hz)]),'ro',label='FFT')
plt.xlabel('频率（奈奎斯特频率的分数）')
plt.ylabel('幅值（任意单位）')
plt.title("手动傅里叶变换以及内置 FFT 函数傅里叶变换的结果对比")
plt.legend()
plt.show()

这段代码展示了如何对一个时间序列信号 x 执行傅里叶变换，以计算其频率成分。代码首先创建一个时间向量 t，随后初始化傅里叶系数数组 fc。对于信号中的每个频率点 f，代码生成一个复数正弦波 csw，并计算该波与信号 x 的点积，从而得到对应的傅里叶系数。同时，使用 NumPy 的快速傅里叶变换 np.fft.fft 方法计算信号的傅里叶变换 fc2，以便进行比较。最后，代码绘制了两种方法得到的傅里叶变换的幅值（避免绘制负频率）

2.3 基于傅里叶变换的步态信息分析

2.3.1 实验数据下载
从网站 https://physionet.org/physiobank/database/gaitdb/下载帕金森病人和对照组的步态信息的数据， 其中帕金森病人的数据文档名称为 pd1-si.txt， 对照组数据文档的名称为 y1-23.si.txt，这里是运用了处理后的数据文件gait.mat。
2.3.2 绘制曲线
通过 Python 编程， 读取两组数据， 绘制帕金森病人和对照组的时域和频域的步态信息对比图，通过观察对比总结从时域和频域步态信息中对帕金森病人进行确诊的标准。

import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack
import scipy.io as sio
import matplotlibmatplotlib.rc("font",family='KaiTi')
matplotlib.rcParams['axes.unicode_minus']=Falsematdata = sio.loadmat('/content/gait.mat')# 从 mat 文件中提取数据
park = matdata['park']
cont = matdata['cont']# 绘制帕金森病人和对照组的步态速度
plt.subplot2grid((3,1),(0,0))
plt.plot(park[:,0],park[:,1],label='帕金森病人')
plt.plot(cont[:,0],cont[:,1],label='对照组')
plt.xlabel('时间（秒）')
plt.ylabel('跨步时间（秒）')
plt.legend()# 定义采样率
srate  = 1000# 创建帕金森病人的步伐时间序列
parkts = np.zeros(int(park[-1,0]*1000))
for i in range(0,len(park)):parkts[int(park[i,0]*1000-1)] = 1# 时间向量和时间点数
parktx = np.arange(0,len(parkts))/srate
parkn  = len(parktx)# 为对照组数据重复上述步骤
contts = np.zeros(int(cont[-1,0]*1000))
for i in range(0,len(cont)):contts[int(cont[i,0]*1000-1)] = 1# 时间向量和时间点数
conttx = np.arange(0,len(contts))/srate
contn  = len(conttx)# 绘制步伐的时间过程
plt.subplot2grid((3,1),(1,0))
plt.plot(parktx,parkts)
#plt.plot(conttx,contts)
plt.xlim([0,50])# 计算两个数据集的功率谱
parkPow = 2*np.abs(scipy.fftpack.fft(parkts)/parkn)
contPow = 2*np.abs(scipy.fftpack.fft(contts)/contn)# 为每个受试者计算单独的频率向量
parkHz = np.linspace(0,srate/2,int(np.floor(parkn/2)+1))
contHz = np.linspace(0,srate/2,int(np.floor(contn/2)+1))# 显示功率谱
plt.subplot2grid((3,1),(2,0))
plt.plot(parkHz[1:],parkPow[1:len(parkHz)])
plt.plot(contHz[1:],contPow[1:len(contHz)])
plt.xlim([0,7])
plt.show()

这段Python代码是为了分析帕金森病患者与健康对照组在步态特征上的差异，并使用频率分析方法来进一步深入理解这些差异。首先，代码通过scipy.io.loadmat函数从一个.mat文件中加载包含帕金森病患者和健康对照组的步态数据。数据被存储在两个不同的变量park和cont中，分别代表帕金森病患者组和对照组的步态信息。
接着，代码在第一个子图中绘制了两组的跨步时间随时间变化的数据，这个对比图直观地展示了两组在步态节奏和速度上的差异。为了进一步的数字信号分析，代码为每个步态事件创建了一个二进制时间序列，在这个序列中，步态事件发生的时间点被标记为1，其他时间则为0。这样的处理将连续的时间事件转换为离散的时间点，便于进行后续的傅里叶变换。
使用傅里叶变换，代码将这些二进制时间序列从时域转换到频域，计算得到两组的功率谱。功率谱显示了不同频率下步态事件的能量密度，通过比较帕金森病患者和健康对照组的功率谱，可以观察到两组在步态控制和协调能力上的明显差异。特别是在较低的频率范围内，这通常与步行节奏和稳定性有关。

2.4 基于傅里叶变换的卷积分析

通过 Python 编程， 绘制激励信号和冲激响应图像。并对激励信号和冲激响应分别进行傅里叶变换，得到它们的频域表示，在频域中将这两个表示相乘得到输出信号的频域表示， 最后通过逆傅里叶变换得到时域的输出信号，同时， 通过时域卷积，得到输出信号，对比两种计算过程得到的结果。

# 信号和卷积核的长度
m = 50  # 信号长度
n = 11  # 卷积核长度# 生成激励信号 (例如，方波信号)
signal = np.zeros(m)
signal[20:30] = 1# 生成分段线性的冲激响应
kernel = np.zeros(n)
x = np.linspace(0, 10, n)
for i in range(n):if 0 <= x[i] < 1:kernel[i] = 0elif 1 <= x[i] <= 2:kernel[i] = x[i] - 1elif 2 < x[i] <= 7:kernel[i] = -0.2 * x[i] + 1.4elif 7 < x[i] <= 10:kernel[i] = 0# 画出信号和冲激响应
plt.figure(figsize=(12, 10))# 绘制激励信号
plt.subplot(4, 1, 1)
plt.plot(signal, label='激励信号')
plt.title('激励信号')
plt.legend()# 绘制冲激响应
plt.subplot(4, 1, 2)
plt.plot(x, kernel, label='冲激响应')
plt.title('冲激响应')
plt.legend()# 对信号和冲激响应进行傅里叶变换
signal_fft = np.fft.fft(signal, m + n - 1)
kernel_fft = np.fft.fft(kernel, m + n - 1)# 在频域中将信号和冲激响应的频域表示相乘
output_fft = signal_fft * kernel_fft# 通过逆傅里叶变换得到时域的输出信号
output_signal_fft = np.fft.ifft(output_fft).real# 直接进行卷积得到的输出信号
output_signal_direct = np.convolve(signal, kernel)# 绘制卷积结果对比
plt.subplot(4, 1, 3)
plt.plot(output_signal_fft, label='频域卷积 + IFFT')
plt.plot(output_signal_direct, 'o', label='时域卷积')
plt.title('卷积结果')
plt.legend()plt.tight_layout()
plt.show()

这段代码主要展示了如何使用卷积来处理信号，通过时域和频域两种不同的方法来实现。首先，代码生成了一个方波作为激励信号和一个分段线性函数作为冲激响应。接着，它绘制了这两个信号的图形。然后，代码通过计算这两个信号的傅里叶变换，并在频域中相乘，最后通过逆傅里叶变换将结果转换回时域，从而得到第一种卷积结果。同时，代码也直接在时域中使用卷积函数进行卷积，以获得第二种卷积结果。这两种结果被绘制在同一图中进行比较。通过这样的处理，代码展示了时域卷积和利用傅里叶变换进行的频域卷积之间的等效性，并可视化了两种方法的卷积输出，验证了频域卷积的正确性和效率。

3.数据处理

3.1数据预处理和参数设置：
①设置中文字体为楷体以确保图表中的中文显示正确，同时设置axes.unicode_minus为False以保证负号显示正确。
②定义了信号的长度为50（变量m）和卷积核的长度为11（变量n），这些参数决定了信号和冲激响应的数据结构。
3.2生成激励信号：
①创建一个长度为50的全零数组，然后在数组的20到30的位置填充1，生成一个方波信号。这个方波信号作为系统的输入，用于后续的卷积操作。
3.3生成分段线性的冲激响应：
①同样创建一个长度为11的全零数组，并通过np.linspace生成一个从0到10等间隔的11个点的数组x，用于定义冲激响应的形状。
②根据x的值，设定冲激响应在不同段的线性关系，生成一个分段线性的响应信号。这个响应信号反映了系统对输入信号的反应。
3.4信号的卷积处理：
①对激励信号和冲激响应进行傅里叶变换。使用np.fft.fft对两个信号进行变换，并确保变换后的长度为信号和冲激响应长度之和减一，以匹配它们的卷积结果的长度。
②在频域中将两个信号的傅里叶变换结果相乘，得到卷积的频域表示。
③使用逆傅里叶变换np.fft.ifft将卷积结果从频域转换回时域，并取实部作为最终的卷积输出信号。
3.5结果的可视化：
①使用matplotlib的subplot功能，在一个图形窗口中创建四个子图。分别绘制激励信号、冲激响应、频域卷积结果及其与直接在时域中计算的卷积结果的比较。
②为每个图添加标题和图例，以清晰地标示每个子图表示的内容。

三、实验结果

1.数据展示

图 1通过自回归过程生成信号

图 2手动傅里叶变换以及内置 FFT 函数傅里叶变换的结果对比

图 3时域和频域的帕金森病人和对照组的步态信息对比

![在这里插入图片描述](https://i-blog.csdnimg.cn/direct/0182f45271c5469ca87a446834c96beb.png

图 4激励信号， 冲激响应和频域和时域分别得到的卷积结果对比

2.结果分析

图2展示了实验信号通过手动傅里叶变换（FT）和快速傅里叶变换（FFT）得到的频域表示。可以观察到，两种方法得到的频率成分在主要的峰值上显示出较高的一致性，验证了FFT方法在实际应用中的有效性和精确度。图中的主要峰值集中在频率0.0、0.2和0.4处，显示了信号的主要频率成分。这一结果对理解信号的频域特性提供了重要的视角，并支持了使用FFT作为信号分析的可靠方法。
图3的顶部子图显示了两组信号的时间序列对比。这些信号代表了帕金森病患者与健康对照组的步态数据，反映了步态速度或节奏的差异。通过比较这些信号的时域和频域数据，可以帮助确诊帕金森病 。
图3的中间子图的二进制信号是对步态事件的标记，即步态开始和结束的时间点，这对于分析步态的节奏和持续性至关重要。这种标记方法有助于在时域分析中准确识别重要事件，从而进行更精确的频域分析。
图3的底部子图的频域分析显示了两个信号的频率成分。在实验中，使用了傅里叶变换来从时域信号中提取频率成分，以分析和比较不同样本的步态特征 。频域分析揭示了帕金森病患者可能存在的步态控制和协调能力的缺失，这通常在低频范围内表现得更为明显。

四、讨论

1.结果解释

在本次实验中，我们通过对帕金森病人和健康对照组的步态数据进行详细的频域和时域分析，观察到帕金森病患者的步态特征在频域中表现出较低的频率成分增强，这与步态的不规则性和减缓有关。这一发现与预期的理论分析相符，支持了使用傅里叶变换分析步态数据作为诊断工具的可行性。然而，实验结果也显示出一些与预期不符的差异，特别是在高频成分的分析上。预期中，高频成分应较低，因为帕金森病患者的快速运动能力受限。但数据显示，部分患者的高频成分异常升高，这可能与个体病情的差异或者是测量误差有关。

2.问题探讨

实验中遇到的主要问题是数据噪声较大，这可能由步态测量设备的精度限制或者是患者在实验中的非标准化运动引起。为解决这一问题，我们采用了多次平均和信号滤波的方法来减少噪声干扰。此外，部分数据处理和分析步骤的计算资源需求较高，限制了实验的效率。针对这一点，我们优化了数据处理算法，采用更高效的数学工具库来处理数据，显著提高了处理速度。

3.改进建议

针对本次实验的经验，我们建议未来研究中可以引入更先进的测量设备，以提高数据的准确性和可靠性。此外，考虑到帕金森病患者步态的多样性，增加样本量并采用分层抽样的方法来考察不同类型和阶段的帕金森病患者，可能会获得更具代表性和普遍适用性的结果。为进一步提高分析的准确性和综合性，可以使用机器学习方法来分析步态数据，通过训练算法自动识别和量化疾病影响的步态模式。这不仅可以提高诊断的精确度，还能在实时监测帕金森病患者的病情进展中发挥重要作用。

参考文献

[1] Mittra Y, Rustagi V. Classification of subjects with Parkinson’s disease using gait data analysis[C]//2018 International Conference on Automation and Computational Engineering (ICACE). IEEE, 2018: 84-89.
[2] Trabassi D, Serrao M, Varrecchia T, et al. Machine learning approach to support the detection of Parkinson’s disease in IMU-based gait analysis[J]. Sensors, 2022, 22(10): 3700.
[3] Shin J H, Yu R, Ong J N, et al. Quantitative gait analysis using a pose-estimation algorithm with a single 2D-video of Parkinson’s disease patients[J]. Journal of Parkinson’s disease, 2021, 11(3): 1271-1283.
[4] Islam J, Ahamed I, Gafur E A. Data Validation using GAIT Signals of Patients to Detect Possibility of Parkinson’s Disease[D]. Department of Electronic and Telecommunication Engineering, 2021.
[5] Sousa J, Silva J, Leonardo R, et al. Inertial-based Gait Analysis Applied to Patients with Parkinson Disease[C]//BIOSIGNALS. 2021: 327-334.
[6] Alharthi A S, Casson A J, Ozanyan K B. Gait spatiotemporal signal analysis for Parkinson’s disease detection and severity rating[J]. IEEE Sensors Journal, 2020, 21(2): 1838-1848.

附录

!pip install numpy
!pip install matplotlib
!pip install scipyimport matplotlib.pyplot as plt
import numpy as np
import matplotlib# matplotlib.rc("font",family='/content/Kaiti.ttf')
# matplotlib.rcParams['axes.unicode_minus']=Falsematplotlib.font_manager.fontManager.addfont('/content/drive/MyDrive/Colab Notebooks/Kaiti.ttf')
matplotlib.rc("font",family='Kaiti')
matplotlib.rcParams['axes.unicode_minus']=False"""3-1"""
# temporal weighting vector
w = np.array([-.6,.9]) # -2,-1
k = len(w)N = 200
x = np.zeros(N)for i in range(k,N):x[i] = sum( w*x[i-k:i] ) + np.random.randn()# add sine wave
x += np.sin(np.linspace(0,10*np.pi,N))# plot the signal
plt.plot(x)
plt.xlabel('X 轴')
plt.ylabel('函数值')
plt.title('通过自回归过程生成信号')
plt.show()"""3-2"""
# time vector
t = np.arange(0,N)/N# initialize Fourier coefficients
fc = np.zeros(N,dtype=complex)for f in range(N):# create complex sine wavecsw = np.exp(-1j*2*np.pi*f*t)# dot product with signalfc[f] = np.dot(csw,x)# check against the fft
fc2 = np.fft.fft(x)# normalized frequencies vector (avoid plotting the negative frequencies)
hz = np.linspace(0,1,int(N/2+1))# and plot
plt.plot(hz,np.abs(fc[:len(hz)]),label='手动 FT')
plt.plot(hz,np.abs(fc2[:len(hz)]),'ro',label='FFT')
plt.xlabel('频率（奈奎斯特频率的分数）')
plt.ylabel('幅值（任意单位）')
plt.title("手动傅里叶变换以及内置 FFT 函数傅里叶变换的结果对比")
plt.legend()
plt.show()import numpy as np
import matplotlib.pyplot as plt
import scipy.fftpack
import scipy.io as sio
import matplotlibmatplotlib.rc("font",family='KaiTi')
matplotlib.rcParams['axes.unicode_minus']=Falsematdata = sio.loadmat('/content/drive/MyDrive/Colab Notebooks/gait.mat')# 从 mat 文件中提取数据
park = matdata['park']
cont = matdata['cont']# 绘制帕金森病人和对照组的步态速度
plt.subplot2grid((3,1),(0,0))
plt.plot(park[:,0],park[:,1],label='帕金森病人')
plt.plot(cont[:,0],cont[:,1],label='对照组')
plt.xlabel('时间（秒）')
plt.ylabel('跨步时间（秒）')
plt.legend()# 定义采样率
srate  = 1000# 创建帕金森病人的步伐时间序列
parkts = np.zeros(int(park[-1,0]*1000))
for i in range(0,len(park)):parkts[int(park[i,0]*1000-1)] = 1# 时间向量和时间点数
parktx = np.arange(0,len(parkts))/srate
parkn  = len(parktx)# 为对照组数据重复上述步骤
contts = np.zeros(int(cont[-1,0]*1000))
for i in range(0,len(cont)):contts[int(cont[i,0]*1000-1)] = 1# 时间向量和时间点数
conttx = np.arange(0,len(contts))/srate
contn  = len(conttx)# 绘制步伐的时间过程
plt.subplot2grid((3,1),(1,0))
plt.plot(parktx,parkts)
#plt.plot(conttx,contts)
plt.xlim([0,50])# 计算两个数据集的功率谱
parkPow = 2*np.abs(scipy.fftpack.fft(parkts)/parkn)
contPow = 2*np.abs(scipy.fftpack.fft(contts)/contn)# 为每个受试者计算单独的频率向量
parkHz = np.linspace(0,srate/2,int(np.floor(parkn/2)+1))
contHz = np.linspace(0,srate/2,int(np.floor(contn/2)+1))# 显示功率谱
plt.subplot2grid((3,1),(2,0))
plt.plot(parkHz[1:],parkPow[1:len(parkHz)])
plt.plot(contHz[1:],contPow[1:len(contHz)])
plt.xlim([0,7])
plt.show()# 信号和卷积核的长度
m = 50  # 信号长度
n = 11  # 卷积核长度# 生成激励信号 (例如，方波信号)
signal = np.zeros(m)
signal[20:30] = 1# 生成分段线性的冲激响应
kernel = np.zeros(n)
x = np.linspace(0, 10, n)
for i in range(n):if 0 <= x[i] < 1:kernel[i] = 0elif 1 <= x[i] <= 2:kernel[i] = x[i] - 1elif 2 < x[i] <= 7:kernel[i] = -0.2 * x[i] + 1.4elif 7 < x[i] <= 10:kernel[i] = 0# 画出信号和冲激响应
plt.figure(figsize=(12, 10))# 绘制激励信号
plt.subplot(4, 1, 1)
plt.plot(signal, label='激励信号')
plt.title('激励信号')
plt.legend()# 绘制冲激响应
plt.subplot(4, 1, 2)
plt.plot(x, kernel, label='冲激响应')
plt.title('冲激响应')
plt.legend()# 对信号和冲激响应进行傅里叶变换
signal_fft = np.fft.fft(signal, m + n - 1)
kernel_fft = np.fft.fft(kernel, m + n - 1)# 在频域中将信号和冲激响应的频域表示相乘
output_fft = signal_fft * kernel_fft# 通过逆傅里叶变换得到时域的输出信号
output_signal_fft = np.fft.ifft(output_fft).real# 直接进行卷积得到的输出信号
output_signal_direct = np.convolve(signal, kernel)# 绘制卷积结果对比
plt.subplot(4, 1, 3)
plt.plot(output_signal_fft, label='频域卷积 + IFFT')
plt.plot(output_signal_direct, 'o', label='时域卷积')
plt.title('卷积结果')
plt.legend()
plt.tight_layout()
plt.show()