MATLAB代码纠错与优化实践

table {
border-collapse: collapse;
width: 100%;
margin-bottom: 1rem;
}
th, td {
border: 1px solid #ddd;
padding: 8px;
text-align: left;
}
th {
background-color: #f2f2f2;
}
tr:nth-child(even) {
background-color: #f9f9f9;
}
pre {
background-color: #f8f8f8;
padding: 15px;
border-radius: 4px;
overflow-x: auto;
}

79、以下两段MATLAB代码段落中，每段至少有一个错误（a部分一个，b部分一个）。找出并纠正这两个错误。a. fftEEG = fft(csd(7,:,10),n_conv); fftWave = fft(wavelets(3,:),n_conv); as = ifft(mean(fftEEG.

fftWave),n_conv); b. sine_wave = exp(2

1i

pi

frequency(10).

time); gaus_win = exp(time.^2./(2

(6/(2

pi

frex(10)))^2)); wavelet = sine_wave .* gaus_win;

以下是将给定文本内容调整为

Markdown 格式

的结果：

a部分错误在于

mean(fftEEG.*fftWave)
，

fftEEG
和

fftWave
是向量，使用

mean
会将向量元素求平均，可能不是预期操作，应直接相乘后进行逆傅里叶变换。

修正后代码：


fftEEG = fft(csd(7,:,10),n_conv);
fftWave = fft(wavelets(3,:),n_conv);
as = ifft(fftEEG.*fftWave,n_conv);

b部分错误在于高斯窗函数

gaus_win
，指数部分应该是负的，以形成高斯形状。

修正后代码：


sine_wave = exp(2*1i*pi*frequency(10).*time);
gaus_win = exp(-time.^2./(2*(6/(2*pi*frex(10)))^2));
wavelet = sine_wave .* gaus_win;

80、在以下两组代码中，每组有一行代码正确，一行代码有错误。找出并修正错误。代码如下：hz = linspace(0,nyquist,floor(N/2)+1); hz = linspace(0,nyquist,floor(N/2)-1); n_conv = size(csd,3) + length(wavelet_time) – 1; n_conv = length(size(csd,3))+length(wavelet_time)–1;

第一组中，

hz = linspace(0,nyquist,floor(N/2)+1);
正确，

hz = linspace(0,nyquist,floor(N/2)-1);
错误，

应改为

hz = linspace(0,nyquist,floor(N/2)+1);
。

第二组中，

n_conv = size(csd,3) + length(wavelet_time) – 1;
正确，

n_conv = length(size(csd,3)) + length(wavelet_time) – 1;
错误，

应改为

n_conv = size(csd,3) + length(wavelet_time) - 1;
（注意原代码中减号为中文符号，应改为英文符号）。

81、下面的代码是正确的，但可能无法产生预期的结果。问题是什么，如何修复？代码为：x = 0:100; gaus = exp(-(x.^2)/100);

问题在于生成的高斯函数可能未按预期归一化，且绘图范围可能无法清晰展示高斯函数特征。可进行归一化处理，同时调整绘图范围。修复代码如下：


x = 0:100;
gaus = exp(-(x.^2)/100);
gaus = gaus / max(gaus);
plot(x, gaus);

归一化操作使高斯函数最大值为1，绘图能直观展示其特征。

82、在这个练习中，你将探究 Morlet 小波的周期数参数控制时间和频率精度权衡的原因。生成五个频率为 10 Hz 的小波，将周期数从 2 线性递增到 15。使用 5×2 的子图布局，在左图中展示对正弦波进行衰减的时域高斯函数，在右图中展示小波功率的频域表示。为每个子图添加标题以表明周期数。这个图能帮助你理解时频精度的权衡吗？（如果不能，请绘制一个更好的图！）

要完成此练习，可按以下步骤操作：

生成五个频率为 10 Hz 的小波，周期数从 2 线性递增到 15；

使用 5×2 的子图布局，左图展示时域高斯函数，右图展示小波功率的频域表示；

为每个子图添加标题表明周期数；

观察图形判断是否能理解时频精度的权衡，若不能则绘制更好的图。

由于没有具体代码实现，无法给出代码示例，但思路是明确的。

83、构建一系列有限脉冲响应（FIR）滤波器（带通：8 – 12 Hz），滤波器阶数从8 Hz的1个周期变化到15个周期，保持其他参数不变。展示滤波器内核在时域的图像（可以给每个内核的y轴添加一个小偏移以提高可见性），以及它们的功率谱在一个按频率 – 阶数排列的图像中（需要进行零填充以使功率谱具有可比性）。滤波器的频率特性如何随阶数变化？接下来，创建一个以10 Hz为中心的实值Morlet小波。将每个滤波器应用于该小波，并测量每个滤波结果的经验半高宽（FWHM）。（提示：在abs(hilbert(filtsig))的结果上测量FWHM）。绘制FWHM随滤波器阶数变化的图像。根据这些结果，在什么阶数点上你会基于不同阶数对结果得出不同的结论？

本题需按以下步骤完成：

构建滤波器

：构建一系列带通为8 – 12 Hz的FIR滤波器，阶数从8 Hz的1个周期到15个周期变化，保持其他参数不变。

绘制图像

：

– 绘制滤波器内核在时域的图像，可添加y轴小偏移；

– 绘制功率谱在频率 – 阶数图像中，进行零填充。

观察滤波器频率特性随阶数的变化规律。

创建小波

：创建以10 Hz为中心的实值Morlet小波。

应用滤波器并测量FWHM

：将每个滤波器应用于小波，在

abs(hilbert(filtsig))
结果上测量FWHM。

绘制FWHM图像

：绘制FWHM随滤波器阶数变化的图像。

得出结论

：根据上述结果，判断在什么阶数点会基于不同阶数对结果得出不同结论。但具体的频率特性变化规律、FWHM数值及不同结论的阶数点需通过实际编程和计算得出。

84、在时域中生成一个sinc函数。使用fft函数并绘制其功率谱。可能绘制结果为空。在命令窗口中检查傅里叶系数，分析出现这种情况的原因。查看sinc函数的公式，思考在x = 0时会发生什么。添加一些代码以确保整个信号包含有限值（可利用isfinite函数），然后检查频率分布。

可按以下步骤完成操作：

在时域生成sinc函数；

使用

fft
函数计算其傅里叶变换并绘制功率谱；

若绘制为空，在命令窗口检查傅里叶系数；

查看sinc函数公式，思考

x = 0
时情况；

使用

isfinite
函数添加代码确保信号含有限值；

检查频率分布。

85、创建一个线性啁啾信号并设计一个窄带通滤波器。从经过滤波的啁啾信号的希尔伯特变换中提取时变功率。同时提取滤波器核的功率谱。在滤波器阶数因子从 1 到 8 的范围内进行循环，重复上述过程。然后将所有功率时程曲线绘制在同一张图上，并在另一张单独的图上绘制所有滤波器核的功率谱。

可按照以下步骤实现该需求：

创建线性啁啾信号。

设计窄带通滤波器。

对啁啾信号进行滤波，并对滤波后的信号进行希尔伯特变换，提取时变功率。

提取滤波器核的功率谱。

设定循环，让滤波器阶数因子从 1 到 8 变化，重复步骤 1 – 4。

将所有功率时程曲线绘制在同一张图上。

在另一张单独的图上绘制所有滤波器核的功率谱。

86、已知线性调频功率时间序列呈高斯状。估算每个滤波器阶数因子对应的功率时间序列的半高宽（FWHM）。绘制半高宽随阶数变化的图。根据估算结果和所绘制的图，你对滤波器阶数参数有什么了解？你是否有办法来验证滤波器参数的合理性？

关于滤波器阶数参数，绘制

半高宽随阶数变化的图

后，可了解到滤波器阶数对滤波效果的影响，如：

不同阶数下滤波器

频率特性的变化

对信号处理后结果（如

半高宽

）的影响

验证滤波器参数合理性的方法有：

对比不同阶数下的滤波结果

，观察是否符合预期的频率特性

与理论计算或已知标准结果进行对比

检查滤波后的信号是否保留了期望的特征

且未引入不合理的噪声或失真

87、变量 step3 和 rmsx 之间有区别吗？（rms 是一个计算均方根的函数，所以这里用 rmsx 作为变量名）

有区别。二者虽都与均方根计算相关，但计算方式不同。step3 是分三步计算，rmsx 是一步计算。

88、变量EEG.data是一个三维矩阵（通道×时间×试验）。你的朋友弗朗西斯想通过对时间（即第二维）求平均来计算均方根（RMS）。以下代码行为什么是错误的，如何修正？你对弗朗西斯在类似情况下避免混淆有什么建议？代码为：rmsf = sqrt(mean(squeeze(EEG.data(30,:,:)).^2,2));

代码错误原因在于，弗朗西斯想对时间（第二维）求平均来计算RMS，但代码里

mean
函数的第二个参数设为

2
，这意味着是对第二个维度（时间）之外的维度求平均，与需求不符。

修正方法是把

mean
函数的第二个参数改为

1
，正确代码为：


rmsf = sqrt(mean(squeeze(EEG.data(30,:,:)).^2, 1));

建议弗朗西斯在处理多维矩阵时，仔细明确每个维度代表的含义，写代码前可以先在纸上列出操作步骤和对应的维度，避免维度使用错误。同时，编写代码后进行简单测试，验证结果是否符合预期。

89、Frances 回来了。她只测量一个电极和一次试验，让事情变简单了，但在 MATLAB 编程上仍有问题。她用代码“rmsf = sqrt(mean(data));”计算均方根（RMS），结果是复数。她的编程错误是什么，为什么会产生复数形式的 RMS？

编程错误在于代码未对数据先平方再求平均，而是直接求平均后开方。若数据中存在负数，直接求平均后的值可能为负，对负数开方就会产生复数形式的 RMS。

正确计算 RMS 应先对数据平方，再求平均，最后开方，即


rmsf = sqrt(mean(data.^2));

90、在双 sinc 函数中，找出所有局部最大值，然后从列表中移除振幅大于 1 或小于 0.1 的峰值。接着移除 x = 40 之后的峰值。绘制图形，以确认你的筛选器是否有效。

可按以下步骤操作：

先找出双

sinc
函数所有局部最大值，使用代码

matlab peeks = find(diff(sign(diff(signal))) < 0) + 1

移除振幅大于 1 或小于 0.1 的峰值，可通过判断

signal(peeks)
的值进行筛选

移除

x = 40
之后的峰值，可根据

x(peeks)
的值筛选

最后绘制图形，以确认筛选器是否有效

91、变量hello是一个5×1的向量。以下代码行输出的大小分别是多少？先思考答案，然后在MATLAB中验证。mean(hello) mean(hello,1) mean(hello,2)

下面是给定的【文本内容】：

mean(hello)
和

mean(hello,1)
的输出大小为 1×1，因为它们计算的是向量的均值，结果是一个标量；

mean(hello,2)
的输出大小为 5×1，因为它按列计算均值，对于 5×1 的向量，每一行有一个均值，共 5 个，形成 5×1 的向量。

92、计算去趋势波动分析（DFA）或去趋势移动平均法（DMA）的初始步骤是对信号进行积分以创建无界形式。如何从无界形式回到原始的有界形式？（提示：积分的逆运算是什么，如何在MATLAB中实现？）

积分的逆运算是微分。在 MATLAB 中，可以使用

diff
函数进行数值微分。例如，若有一个向量

y
表示无界形式的信号，

dy = diff(y)
可得到其微分结果，以此尝试从无界形式回到有界形式。

但需注意

diff
函数会使向量长度减 1，可能需要根据具体情况处理边界问题。

93、以下两行代码有什么区别？是否存在区别？何时会有区别？dt = mean(diff(timevec)) / 1000; dt = mean(timevec(2)-timevec(1)) / 1000;


第一行代码中，`diff(timevec)` 计算 `timevec` 相邻元素的差值，`mean(diff(timevec))` 求这些差值的平均值，最终 `dt` 是这个平均值除以 1000。  

第二行代码中，`timevec(2)-timevec(1)` 只是 `timevec` 第二个元素与第一个元素的差值，`mean` 对这个单一差值求平均（结果还是该差值），最终 `dt` 是这个差值除以 1000。  

当 `timevec` 中相邻元素差值都相等时，两行代码结果相同；当 `timevec` 中相邻元素差值不相等时，两行代码结果不同。

94、这里有一种计算尖峰 – 场相干性的新方法：将局部场电位（LFP）相角时间序列离散成 N 个区间（例如 10 个区间），并计算每个区间的平均尖峰数。然后将结果绘制成以尖峰计数为 y 轴、LFP 相位为 x 轴的条形图。对几个频率进行此操作。请简要概括该方法的步骤。

该方法步骤为：

将 LFP 相角时间序列离散成若干区间；

计算每个区间的平均尖峰数；

绘制以尖峰计数为 y 轴、 LFP 相位为 x 轴的条形图；

对几个频率重复上述操作。

95、在计算尖峰节律性时，不使用代码行 spRhy = spRhy([end – win+1:end 1:win+1])，而是使用函数 fftshift 来获得相同的结果，应该怎么做？

可以使用如下代码替代原代码以获得相同结果：


spRhy = fftshift(spRhy);

fftshift
函数会将零频分量移到频谱中心，对于周期性数据，能实现与原代码类似的索引重排效果。

96、路易斯试图细致地检查叠加在同一图上的单试验数据，但MATLAB渲染图形花费了很长时间，部分原因是他使用的是一台旧笔记本电脑（2014年的）。他向你寻求帮助。根据他下面的代码，你至少可以给他哪两条不同的建议？（变量alldata包含60000次试验和1200个时间点）for i=1:size(alldata,1) figure(1), hold on plot(i+alldata(i,:)) end

建议一：不要逐行运行代码，可考虑批量处理数据后再绘图。

建议二：减少绘图次数，比如对数据进行抽样，减少参与绘图的试验数量。

97、有代码行 tempsps(tempsps < – win | tempsps>win) = []; ，不使用“或”运算符重写这一行。（提示：如何使负值和正值相等？）

可以将代码改写为


tempsps(abs(tempsps) > win) = [];

利用绝对值函数

abs()
去掉了“或”运算符。

98、以下代码行以错误的方式减去了尖峰时间序列的均值。你如何判断呢？（提示：合理性检查！）spikes = bsxfun(@minus,spikes,mean(spikes,1));

可以通过合理性检查判断。一般来说，减去均值是为了让数据中心化，即数据围绕零值分布。若减去均值的方式错误，可能导致数据没有正确中心化。

在代码中，

mean(spikes,1)
计算的是每列的均值。若尖峰时间序列应按行计算均值，这样就会出错。

可通过检查减去均值后的数据分布是否合理，比如数据是否大致围绕零值波动，来判断代码是否正确。

99、调整 k – 均值聚类以提取并绘制三个聚类。你是否认同这些结果？这对于 k – 均值聚类和预先知晓聚类数量有什么启示？

可以使用MATLAB的

kmeans
函数进行调整，示例代码如下：


[gidx, cents, sumdist, distances] = kmeans(d, 3);

其中

d
是 75 行 2 列的矩阵，包含 75 个数据点的 x 和 y 值，输入

3
指示提取三个聚类。

之后可通过绘图代码直观展示聚类结果。是否认同结果需结合具体数据情况判断。若结果符合数据的实际分布特征，则可认同。

这表明 k-均值聚类结果依赖于预先设定的聚类数量，若能事先知晓合适的聚类数量，聚类结果会更可靠；若缺乏先验知识，可能导致聚类结果不稳定或不合理。

100、现在你有两个尖峰簇，测试这些单元是否相互作用。识别簇2中在簇1中的尖峰前3毫秒或后3毫秒产生尖峰的动作电位。绘制它们的波形形状。为作比较，同时绘制随机选择的尖峰的平均波形形状。

要完成该任务，需按照以下步骤进行：

筛选符合条件的动作电位

从数据中找出

簇2

中在

簇1

尖峰前后

3毫秒

内产生尖峰的动作电位。

绘制波形

对筛选出的动作电位绘制其波形图。

随机选择与计算平均波形

随机选择部分尖峰，计算其平均波形，并进行绘制。

可使用的工具包括：

–

Python

：NumPy、Matplotlib 库

–

MATLAB

101、已知协方差矩阵（设为变量spikecov）是从数据矩阵spikes计算得到，计算方式为spikecov = spikes’ * spikes ，协方差矩阵spikecov的秩为11，其维度是12×12。那么spikes的秩是多少？基于这个结果，你能推测出关于矩阵A的秩和矩阵ATA的秩之间关系的线性代数法则吗？

由于协方差矩阵

spikecov
是由

spikes
矩阵通过

spikecov=spikes′⋅spikesspikecov=spikes′⋅spikes

计算得到，即 $ A^T A $ 的形式。

已知

spikecov
的秩为 11，维度为 $ 12 imes 12 $，根据线性代数中

rank(ATA)=rank(A)rank(ATA)=rank(A)

的法则，可推测

spikes
矩阵的秩为 11。

线性代数法则为

：

矩阵 $ A $ 的秩等于矩阵 $ A^T A $ 的秩，即

rank(ATA)=rank(A)rank(ATA)=rank(A)

102、编写额外的查看代码，以 3×4 的子图矩阵显示结构 MRI 扫描在第一个和最后一个切片之间的 12 个线性间隔切片。生成对应三个维度的三幅图。

以下是实现该功能的 MATLAB 代码示例：


% 读取结构 MRI 数据
strMRI = readnifti('MNI152_T1_1mm.nii');

% 获取每个维度的大小
[sizeX, sizeY, sizeZ] = size(strMRI);

% 生成 12 个线性间隔的切片索引
sliceIndices = linspace(1, sizeX, 12);

% 第一个维度的图
figure;
for i = 1:12
    sliceIndex = round(sliceIndices(i));
    subplot(3, 4, i);
    imagesc(squeeze(strMRI(sliceIndex, :, :))');
    axis image;
    axis xy;
    colormap gray;
end

% 生成 12 个线性间隔的切片索引
sliceIndices = linspace(1, sizeY, 12);

% 第二个维度的图
figure;
for i = 1:12
    sliceIndex = round(sliceIndices(i));
    subplot(3, 4, i);
    imagesc(squeeze(strMRI(:, sliceIndex, :))');
    axis image;
    axis xy;
    colormap gray;
end

% 生成 12 个线性间隔的切片索引
sliceIndices = linspace(1, sizeZ, 12);

% 第三个维度的图
figure;
for i = 1:12
    sliceIndex = round(sliceIndices(i));
    subplot(3, 4, i);
    imagesc(squeeze(strMRI(:, :, sliceIndex))');
    axis image;
    axis xy;
    colormap gray;
end

这段代码首先读取结构 MRI 数据，然后分别针对三个维度，生成 12 个线性间隔的切片索引，在每个维度上创建一个 3×4 的子图矩阵来显示这些切片。

103、结构扫描具有高分辨率。尝试通过创建包含跳过体素的新图像来对图像进行下采样（例如，每隔一个体素取一个，或每隔两个体素取一个，以此类推）。计算每个下采样版本的大小相对于原始版本大小的百分比。在不同下采样图像的z平面上绘制最接近50%的切片，并在图像标题中注明尺寸减小的比例。

本题是一个图像处理任务，可按以下步骤完成：

对图像进行下采样，创建包含跳过体素的新图像；

计算每个下采样版本的大小相对于原始版本大小的百分比；

在不同下采样图像的 z 平面上找到最接近 50% 的切片；

绘制这些切片，并在图像标题中注明尺寸减小的比例。

104、利用功能磁共振成像（fMRI）数据，计算每个体素随时间变化的信号强度值的均值和标准差。然后在一个图中展示这两个图的一些切片。若想获得额外加分，不使用均值（mean）和标准差（std）函数（或方差 var 函数，用它就太简单了）完成此练习。

可按以下步骤操作：

读取 fMRI 数据；

不使用

mean
和

std
函数计算每个体素信号强度值随时间的均值和标准差；

从均值和标准差图中选取一些切片；

将这些切片展示在一个图中。

105、已知阈值图颜色与代码所写相反（黑色为1，白色为0），通过修改代码行contourf(threshmap)重现此反转图（不要简单地重新定义颜色轴）。有一种添加一个字符的方法和一种添加两个字符的方法，请问分别是什么？

添加一个字符的方法可能是在

threshmap
前加一个波浪号

"~"
，即

contourf(~threshmap)
；添加两个字符的方法可能是在

threshmap
前加

"1 - "
，即

contourf(1 - threshmap)
。

106、下面这行代码不会产生 MATLAB 错误，但可能无法产生预期结果。请说明错误是什么、有什么影响以及如何修复。（dbmap 是一个二维矩阵，代表经过预刺激基线分贝标准化后的时频功率）imthresh = abs(dbmap>3);

错误：

dbmap>3
会生成一个逻辑矩阵，其中大于 3 的元素为

true
（即 1），小于等于 3 的元素为

false
（即 0），对逻辑矩阵取绝对值没有实际意义，无法得到预期结果。

影响：不能对

dbmap
中大于 3 的元素进行正确处理，只是简单得到一个逻辑矩阵的绝对值，值仍为 0 和 1。

修复方法：若要获取

dbmap
中大于 3 的元素的绝对值，可使用


imthresh = abs(dbmap) .* (dbmap>3);

107、创建一个正弦波，使用bwconncomp和bwlabel来确定正弦波中值大于0.5的位置。然后在正弦波中所有值大于0.5的时间点上绘制星号。最后，调整代码，使星号绘制在y = 0直线上。

代码实现正弦波处理

本题要求编写代码完成一系列操作。首先创建正弦波，然后用

bwconncomp
和

bwlabel
函数确定正弦波中值大于0.5的位置，接着在这些位置的正弦波上绘制星号，最后将星号绘制在y = 0直线上。可按以下步骤编写代码：

创建正弦波

找出正弦波中值大于0.5的位置

使用

bwconncomp
和

bwlabel
函数

在正弦波上绘制星号

调整代码将星号绘制在y = 0直线上

示例代码


% 创建正弦波
t = 0:0.01:10;
y = sin(t);

% 找出正弦波中值大于0.5的位置
idx = y > 0.5;

% 使用bwconncomp和bwlabel函数
CC = bwconncomp(idx);
labels = bwlabel(idx);

% 在正弦波上绘制星号
figure;
plot(t, y);
hold on;
plot(t(idx), y(idx), '*');
hold off;

% 调整代码将星号绘制在y = 0直线上
figure;
plot(t, y);
hold on;
plot(t(idx), zeros(size(t(idx))), '*');
hold off;