MATLAB自定义函数

简介

本文总结了一部分我在科研工作中常用的一些自定义的 MATLAB 函数，每个函数均有详细的帮助文档和注释，可供大家参考。详细内容可参考右侧的目录，同时每个函数均提供源代码以及相应的 .m 文件。将该 .m 文件放到你的 MATLAB 的工作目录下即可直接调用之，而且输入 help 函数名 命令还可以查看对应函数的帮助文档。

函数列表

fft_plot —— 快速绘制信号的频谱图

• 点击下载 .m 文件
• 该函数首先对输入信号进行傅里叶变换，然后计算其单边频谱图，随后自动调用 plot 函数进行绘图。同时可选择是否在信号结尾补零（以提高频谱分辨率，参考 快速傅里叶变换(FFT)中为什么要“补零”？ - 知乎)，以及是否在做傅里叶变换变换前对信号加窗（以抑制频谱泄漏和旁瓣，参考 Matlab FFT 加窗 - 知乎），默认选项是补零至 2 的整数次幂，以及使用汉明窗。源代码参考如下：

function [f, P1] = fft_plot(X, Fs, varargin)
% FFT_PLOT 自动使用 fft 计算频谱并随即调用 plot 绘制频谱图
%
% 使用方法
%   fft_plot(X, Fs) 绘制采样率为 Fs 的信号 X 的频谱图
%   fft_plot(X, Fs, 'is_hamming', 1) 是否加窗，默认为 1（加窗），0 则不加
%   fft_plot(X, Fs, 'is_padding', 1) 是否补零，默认为 1（补零），0 则不补
%   [f, P1] = fft_plot(X, Fs) 返回绘图所用的频率 f 和单边频谱 P1，供自定义绘图
%
% 输入参数
%   X  输入信号
%   Fs 采样率
%   is_hamming （可选）是否加窗，默认加窗
%   is_padding （可选）是否补零，默认补零
%
% 输出参数（可忽略）
%   f  绘图所使用的横轴频率
%   P1 绘图所使用的纵轴单边频谱

    p = inputParser;                  % 解析参数
    addRequired(p, 'X');              % 必有参数 X
    addRequired(p, 'Fs');             % 必有参数 Fs
    addOptional(p, 'is_hamming', 1);  % 默认加窗
    addOptional(p, 'is_padding', 1);  % 默认补零
    parse(p, X, Fs, varargin{:});

    w = ones(size(p.Results.X));      % 默认矩形窗
    L = length(p.Results.X);          % 信号的长度
    if p.Results.is_hamming           % 是否加窗
        w = reshape(hamming(L), size(p.Results.X));
    end
    if p.Results.is_padding           % 是否补零
        L = 2 ^ (nextpow2(L) + 1);
    end

    Y = fft(w .* p.Results.X, L);     % 信号的傅里叶变换
    P2 = abs(Y/L);                    % 双侧频谱
    P1 = P2(1:L/2+1);                 % 单侧频谱
    P1(2:end-1) = P1(2:end-1) * 2;
    f = p.Results.Fs * (0:L/2)/L;     % 频率 f
    plot(f, P1);                               % 按照给定参数绘图
    title('Single-Sided Amplitude Spectrum');  % 默认标题
    xlabel('f (Hz)');                          % 默认 xlabel
    ylabel('|P1(f)|');                         % 默认 ylabel

end

light_spot —— 提取图像的荧光点信号

• 点击下载 .m 文件
• 该函数主要通过形态学开运算来提取荧光成像图中的荧光亮点信号，同时调整参数之后也可以用来提取图像中较为明亮的前景信号以去除背景噪声。源代码参考如下：

function img_filtered = light_spot(img, disk_radius, threshold)
% LIGHT_SPOT 提取荧光图像中的光点
%
% 使用方法
%   img_filtered = light_spot(img, disk_radius, threshold)
%
% 输入参数
%   img         输入图像
%   disk_radius 执行形态学开运算的结构单元半径，其值越小。保留结构越多
%   threshold   图像增强卡的阈值，按百分比计算，小于该值的元素将被设置为 0
%
% 输出参数
%   img_filtered  图像增强后的图像

    img2 = imtophat(img, strel('disk', disk_radius));
    img3 = imadjust(img2, [0, 0.1], [0, 1]);
    img_filtered = img3 .* uint8(img3 > 255 * threshold);

end

generate_plane_signal —— 产生无限大平板的光声仿真信号

• 点击下载 .m 文件
• 该函数产生 $\delta$ 脉冲激发的无限大平板的光声信号，使用的仿真公式如下所示：

$p(z,t)=\frac{p_0}{2}\left[u(z-ct+\frac{d}{2})u(-z+ct+\frac{d}{2})+u(z+ct+\frac{d}{2})u(-z-ct+\frac{d}{2})\right]$

function signal = generate_plane_signal(D, r, vs, t)
% GENERATE_PLANE_SIGNAL 产生无线大平板的光声信号
% 其产生的信号的默认初始声压 p0 = 1
%
% 使用方法
%   signal = generate_plane_signal(D, r, vs, t)
%
% 输入参数
%   D  平板的厚度
%   r  探测器距离平板中心的距离
%   vs 声速
%   t  时间序列
%
% 输出参数
%   signal 产生的信号

    u = @(x) double(x > 0);  % u 为单位阶跃函数
    signal = (u(r - vs * t + D / 2) .* u(-r + vs * t + D / 2) + ...
        u(r + vs * t + D / 2) .* u(-r - vs * t + D / 2)) / 2;

end

generate_sphere_signal —— 产生球的光声仿真信号

• 点击下载 .m 文件
• 该函数产生 $\delta$ 脉冲激发的球的光声信号，使用的仿真公式如下所示：

$p(r,t)=p_0\frac{r-ct}{2r}u(R-|r-ct|)$

function signal = generate_sphere_signal(R, r, vs, t)
% GENERATE_SPHERE_SIGNAL 产生球的光声信号
% 其产生的信号的默认初始声压 p0 = 1
%
% 使用方法
%   signal = generate_sphere_signal(R, r, vs, t)
%
% 输入参数
%   R  球半径
%   r  探测器距离球心的距离
%   vs 声速
%   t  时间序列
%
% 输出参数
%   signal 产生的信号

    u = @(x) double(x > 0);  % u 为单位阶跃函数
    signal = (r - vs * t) / (2 * r) .* u(R- abs(r - vs * t));

end

generate_sphere_signal_diff —— 产生球的光声仿真信号的时间导数

• 点击下载 .m 文件
• 该函数产生 $\delta$ 脉冲激发的球的光声信号的时间导数。源代码参考如下：

function signal_diff = generate_sphere_signal_diff(R, r, vs, t)
% GENERATE_SPHERE_SIGNAL 产生球的光声信号的时间导数信号
% 其产生的信号的默认初始声压 p0 = 1
%
% 使用方法
%   signal_diff = generate_sphere_signal(R, r, vs, t)
%
% 输入参数
%   R  球半径
%   r  探测器距离球心的距离
%   vs 声速
%   t  时间序列
%
% 输出参数
%   signal_diff 产生的时间导数信号

    u = @(x) double(x > 0);  % u 为单位阶跃函数
    signal_diff = -vs / (2 * r) .* u(R- abs(r - vs * t));

end

generate_cylinder_signal —— 产生无限长圆柱的光声仿真信号

• 点击下载 .m 文件
• 该函数产生 $\delta$ 脉冲激发的无限长圆柱的光声信号。源代码参考如下：

function signal = generate_cylinder_signal(R, r, vs, t)
% GENERATE_CYLINDER_SIGNAL 产生无限长圆柱的光声信号
% 其产生的信号的默认初始声压 p0 = 1
%
% 使用方法
%   signal = generate_cylinder_signal(R, r, vs, t)
%
% 输入参数
%   R  圆柱的半径
%   r  探测器距离圆柱中心的距离
%   vs 声速
%   t  时间序列
%
% 输出参数
%   signal 产生的信号

    N = length(t);
    signal = zeros(1, N);
    for i = 1:N
        ti = t(i);
        d = vs * ti;
        if ti <= (r - R) / vs
            signal(i) = 0;
        elseif ti > (r - R) / vs && ti < (r + R) / vs
            low = (r - R) / d;
            high = (r^2 + d^2 - R^2) / (2 * r * d);
            y = @(x) atan(abs(sqrt((R^2 - (r - d * x).^2) ./ ...
                (r^2 + d^2 - 2 * r * d * x - R^2))));
            signal(i) = 4 * d^2 * integral(y, low, high) + ...
                pi * d / r * (R^2 - (r - d)^2);
        elseif ti >= (r + R) / vs
            low = (r - R) / d;
            high = (r + R) / d;
            y = @(x) atan(abs(sqrt((R^2 - (r - d * x).^2) ./ ...
                (r^2 + d^2 - 2 * r * d * x - R^2))));
            signal(i) = 4 * d^2 * integral(y, low, high);
        end
        signal(i) = signal(i) / d;
    end
    temp = diff(signal) / (t(2) - t(1)) / (4 * pi * vs^2);
    signal(1:N - 1) = temp;
    signal(N) = signal(N - 1);

end

recon_pact_signal —— 使用 DAS 算法重建光声信号

• 点击下载 .m 文件
• 该函数使用 delay and sum（DAS）算法来重建光声信号，只要按给定的参数格式输入即可。将之微改之后即可用于 GPU 加速的 CUDA 代码的原始代码进行编写，是一个高度结构化的 DAS 算法的实现。源代码参考如下：

function signal_recon = recon_pact_signal( ...
    signal_backproj, detector_location, x_grid, y_grid, z_grid, ...
    vs, fs, num_detector, num_time)
% RECON_PACT_SIGNAL 使用 delay and sum (DAS) 算法重建光声信号
%
% 使用方法
%   signal_recon = recon_pact_signal(signal_backproj, detector_location, ...
%                                    x_grid, y_grid, z_grid, vs, fs, ...
%                                    num_detector, num_time)
% 输入参数
%   signal_backproj         反投影信号，其每一行是一个探测器的信号
%   detector_location       探测器坐标，其每一行是一个探测器的坐标
%   x_grid、y_grid、z_grid  重建区域空间网格
%   vs、fs                  声速、采样率
%   num_detector、num_time  探测器数量、采样点数
%
% 输出参数
%   signal_recon  重建的信号，大小于 x_grid 相同

    signal_recon = 0 * x_grid;  % 初始化重建信号
    for i = 1:num_detector
        temp_signal = signal_backproj(i, :);      % 取出本次循环要用的信号
        temp_location = detector_location(i, :);  % 取出本次循环要用的坐标
        temp_signal(num_time) = 0;                % 方便处理过大的索引
        % 计算距离和延迟索引
        dx = x_grid - temp_location(1);
        dy = y_grid - temp_location(2);
        dz = z_grid - temp_location(3);
        d = sqrt(dx.^2 + dy.^2 + dz.^2);
        idx = round(d / vs * fs);
        % 处理过大的索引
        idx(idx > num_time) = num_time;
        % delay and sum
        if idx < num_time
            signal_recon = signal_recon + temp_signal(idx);
        end
    end

end

To Be Continued…