Matlab编码说明

Question

我在互联网上找到了这个 Matlab 代码。它实际上是整个代码的一部分，可以在here 找到。有人可以逐行解释发生了什么。我真的很绝望……

% ------------- % This is code to make the edge detecting filter % ----%
function filter=gaussfilt(N)

% calculate alpha so the filter fills all N points
alpha=N;
first=-(1-N/2)*exp(-(1-N/2)^2/alpha);
count=0;
while first<.1*(-(1530/4000*N-N/2)*exp(-(1530/4000*N-N/2)^2/alpha))
    count=count+1;
    alpha=N*500*count;
    first=-(1-N/2)*exp(-(1-N/2)^2/alpha);
end

for n=1:N
     filter(n)=-(n-N/2)*exp(-(n-N/2)^2/alpha);   % d/dt of a gaussian
end
filter=filter/sum(abs(filter));     % normalization

return

Answer 1

此函数尝试返回高斯梯度 - 当您将其与一组 o 数据点进行卷积时，它将检测数据中的边缘，同时对平均水平没有突然变化的区域中的点进行平滑处理。这是 N = 15 时的输出：

但是，如果将 N 增加到 20，就会发生一些疯狂的事情，因为代码坦率地说是错误的（也很丑陋）。你最终得到的不是一条漂亮的曲线，而是一条直线——一个糟糕的过滤器。这是因为为 alpha 计算新值的尝试非常糟糕。

Robert P 已经提供了代码正在执行的逐步描述。让我向您展示编写此函数的“正确方法”（您会看到我使用了 Robert 提到的一些技术，另外其他一些我会解释的）...

function myFilter = gaussfilt2(N, alpha)
% myFilter = gaussfilt(N, alpha)
% returns an N point normalized array of filter coefficients
% corresponding to the gradient of a Gaussian over the interval [-1 1]
% with a standard deviation of 1/alpha 
% in other words, the higher alpha, the sharper the filter
% default value for alpha is 3

if ~exist('alpha', 'var')
  alpha = 3; 
end 

x = linspace( -1, 1, N); % create a vector of N values between -1 and 1 inclusive
sigma = 1.0 / alpha; % convert from alpha to sigma as used in Gaussian formula

% compute first derivative, but leave constants out
% we will normalize later by summing over the coefficients
myFilter = -x .* exp( -(x.^2)/(2*sigma.^2)); % using .* for element-by-element operation

% normalize:
myFilter = myFilter / sum( abs( myFilter ) ); % absolute sum of coefficients is now one

当您使用N=150 和alpha = 3 运行此函数时，曲线如下所示：

一些需要解释的技术：

文档当您在函数声明正下方开始注释块（用% 表示）时，该块中的所有内容都会在用户请求此函数的帮助时显示给用户（例如，通过突出显示另一个脚本中的函数名称并按 F1)。总是个好主意

矢量化 Matlab 不擅长做“显式循环”，但非常擅长“隐式循环”。只要有可能，您就想一次对一大堆值进行“相同的计算”。在这种情况下，我计算一个向量 xusing thelinspace( x1, x2, n )function. This returns an equally spaced array ofnvalues fromx1tox2` 包括在内。现在我可以用一条语句计算整个函数

逐元素乘法 Matlab 真正用于矩阵操作（即“Matlab”中的“Mat”——它代表矩阵，而不是数学）。如果您有两个向量 a 和 b，并且您希望将它们逐个元素相乘（因此结果为 [a(1)8b(1) a(2)*b(2) ... a(n)*b(n)]），则使用 .* 运算符。

灵活性 不是“硬连线”alpha 的值，而是将其作为第二个参数允许您重复使用相同的函数并更改滤镜的锐度。允许用户省略变量并提供默认值意味着当函数只有一个参数时编写的程序仍然可以使用它。当名称为 'alpha' 的变量 ('var') 存在时，exist('alpha', 'var') 调用将返回 true。在函数前面添加 ~ 会否定结果 - 就像在某些其他语言中的 If Not。

Answer 2

我完全同意 Floris 的观点，但我会给你一些关于正在发生的事情的线索。根据您的问题，我假设您以前从未使用过 MATLAB，因此我强烈建议您尝试一些在线提供的教程。

function filter=gaussfilt(N)

这定义了函数的名称gaussfilt、输入变量N 和输出变量filter。您的 m 文件必须保存为 gaussfilt.m。

filter 是 MATLAB 中的许多内置函数之一，因此不是一个好的变量名。我建议您使用诸如gauss_filter 之类的名称，或者基本上其他任何名称。在您的情况下，这可能并不重要，因为您没有使用 filter 函数，但尽管如此，使用这样的名称是一个坏习惯。这也适用于size、length、max 等名称。

alpha=N;
first=-(1-N/2)*exp(-(1-N/2)^2/alpha);
count=0;

这些行只是为变量名赋值。 exp() 是自然指数。 count 用于跟踪while 循环循环的次数。

while first<.1*(-(1530/4000*N-N/2)*exp(-(1530/4000*N-N/2)^2/alpha))
    count=count+1;
    alpha=N*500*count;
    first=-(1-N/2)*exp(-(1-N/2)^2/alpha);
end

只要满足以下条件，就会计算while 和end 之间的所有内容。条件应该是不言自明的。

first < .1*(-(1530/4000*N-N/2)*exp(-(1530/4000*N-N/2)^2/alpha))

count 递增，以便您可以跟踪循环执行了多少次。但是，它从未使用过，因此在这种情况下是不必要的。但是，我建议您保留计数器，并在您的 while 循环中包含另一个条件，即如果 count > 1e6 或其他一些较大的数字，它应该停止运行。这样一来，您将避免在不满足其他条件的情况下永远运行循环。

for n=1:N
     filter(n)=-(n-N/2)*exp(-(n-N/2)^2/alpha);   % d/dt of a gaussian
end

for n = 1:N 是一个循环，将运行N 次，其中n 将是1 第一次，2 第二次等等。filter(n) = .. 将值分配给nth 位置在filter 变量中，从而创建一个长度为 N 的向量。这是创建向量的糟糕方式。您应该始终为向量分配内存，以避免向量在循环内增长。 “增长”的向量非常非常慢。因此，在开始循环之前，您应该执行以下操作：

filter = zeros(1,N);

这会创建一个零向量。这可能有点过分，但为filter 赋值的更好方法是使用arrayfun：

filter = arrayfun(@(n) (-(n-N/2) * exp(-(n-N/2)^2 / alpha)), 1:N);

查看this answer 中的第 6 点以了解原因。

最后一行：

filter=filter/sum(abs(filter));

好吧，您将filter 值除以这些值的总和，从而创建一个新的filter，其总和等于1。这是您将从函数中获得的输出。

在使用函数时，必须写如下：

filter_vector = gaussfilt(N) % where N is an integer

最后，请使用空格！它更容易阅读！

再次，我推荐一些 MATLAB 教程...