Simulink仿真运行时输入信号频率估计技术问题咨询

背景描述 #

需要在Simulink仿真运行时（非Matlab后处理环节）实现输入信号的频率估计，编写了如下Matlab Function代码：

function frq_est = fcn(u)

persistent x_buff

num_data = 2000;

if isempty(x_buff),x_buff=zeros(num_data,1);end

x_buff(2:num_data) = x_buff(1:num_data-1);
x_buff(1)          = u;

Nfft    = 1024;
[Pxx,f] = pwelch(x_buff, gausswin(Nfft), Nfft/2, Nfft,1000);

[~,loc] = max(Pxx);

frq_est = f(loc);

问题现象 #

1. 模型在Matlab R2023a中可运行但存在频率估计误差；在R2020a中完全无法运行，报错信息：

Data 'frq_est' is inferred as a variable size matrix, while its properties in the Model Explorer specify its size as inherited or fixed. Please check the 'Variable Size' check box and specify the upper bounds in the size field.
Component:MATLAB Function | Category:Coder error
Simulink cannot determine sizes and/or types of the outputs for block 'frq_est_rt/MATLAB Function' due to errors in the block body, or limitations of the underlying analysis. The errors might be inaccurate. Fix the indicated errors, or explicitly specify sizes and/or types for all block outputs.

2. 频率估计存在恒定偏移：输入20Hz正弦波时估计值为19.5Hz；输入10Hz正弦波时为9.5Hz。

咨询问题（按优先级） #

1. 是否存在工业界广泛使用的、更成熟可靠的Simulink仿真运行时信号频率估计方案？
2. 如何修改现有代码使其在Matlab R2020a中正常运行？
3. 频率估计结果出现恒定偏移的原因是什么？

1. 工业界常用的Simulink运行时频率估计方案 #

工业界成熟的方案主要分为三类：

• 过零检测法：针对正弦类信号，通过检测相邻过零点的时间差计算频率，用Simulink基础模块（比较器、定时器、触发子系统）即可搭建，无需编写代码，实时性强，适合低资源场景。
• 锁相环(PLL)方案：Simulink的Control System Toolbox或Signal Processing Toolbox提供现成的Phase-Locked Loop模块，能稳定跟踪时变频率，抗干扰能力强，适配工业噪声环境。
• 固定点数FFT模块方案：若坚持用频谱分析，可直接使用Simulink的FFT固定点数模块替代Matlab Function中的pwelch，避免代码生成的尺寸推断问题，同时可配合峰值插值提升精度。

2. 适配Matlab R2020a的代码修改 #

R2020a的代码生成对变量尺寸/类型的推断更严格，需明确指定所有关键变量的属性，修改后的代码如下：

function frq_est = fcn(u)
%#codegen
% 提前初始化输出为标量，明确类型
frq_est = 0.0;

persistent x_buff
num_data = 2000;

if isempty(x_buff)
    % 明确指定缓冲区的尺寸和数据类型
    x_buff = zeros(num_data, 1, 'double');
end

% 简化缓冲区更新逻辑，避免索引歧义
x_buff(2:end) = x_buff(1:end-1);
x_buff(1) = u;

Nfft = 1024;
% 单独生成窗函数，明确其尺寸
win = gausswin(Nfft);
% 调用pwelch时固定参数，确保输出频率向量尺寸固定
[Pxx,f] = pwelch(x_buff, win, Nfft/2, Nfft, 1000);

% 强制转换索引为整数类型，避免尺寸推断错误
[~,loc] = max(Pxx);
loc = cast(loc, 'uint32');

% 直接索引频率向量，输出标量结果
frq_est = f(loc);

同时需完成Simulink块配置：

• 打开Model Explorer，找到该Matlab Function块的输出frq_est，勾选Variable Size，设置尺寸上限为1（标量）；
• 确认输入u的类型为double标量，与代码定义一致。

3. 频率估计恒定偏移的原因 #

偏移核心来自FFT的频率分辨率限制与频谱泄漏：

• 采样率设置为1000Hz，Nfft=1024，FFT的频率分辨率为Fs/Nfft = 1000/1024 ≈ 0.9766Hz，即FFT的频率点是离散的0, 0.9766, 1.9531,..., 499.0234Hz；
• 输入的10Hz、20Hz信号并不落在这些离散频率点上，max(Pxx)找到的是最接近的频率点，再加上高斯窗主瓣较宽，会进一步模糊峰值位置，导致估计值偏向更低的频率点；
• 你观察到的9.5Hz、19.5Hz偏移，还与缓冲区长度（2000点）和pwelch的重叠设置有关，重叠数据的频谱叠加进一步拉低了峰值对应的频率位置。

内容的提问来源于stack exchange，提问作者Eric