信号增强算法

Question

我需要一种算法（最好使用类似 Pascal 的语言，但它的结尾并不重要），它将使左侧的“信号”（实际上是一系列数据点）看起来像右侧的。

信号来源：
信号由机器产生。解释过于简单，机器正在测量流过透明管的液体的密度。因此，该信号与电信号（音频/射频）完全不同。 数据点可能如下所示： [1, 2, 1, 3, 4, 5, 4, 3, 2, 1, 13, 14, 15, 18, 23, 19, 17, 15, 15, 15, 14 , 11, 9, 4, 1, 1, 2, 2, 1, 2]

我需要什么：
我需要准确检测“峰值”。为此，我已经有一段代码，但它不适用于下图所示的不良信号。
我认为我们可以将其视为意外通过低通滤波器的信号，现在我想恢复它。

注意事项：
有 4 个信号，但它们是分开的，因此可以单独分析。因此，只需考虑如何处理其中一个就足够了。

在一个峰值之后，如果信号下降得不够快，我们可以认为存在多个峰值（您可以在系列末尾的“红色”信号中最好地看到这一点）。

优点是全系列都可用（所以信号不是实时的，已经存档了）！

[由 Spektre 编辑]我从图像中提取了红色样本点

float f0[]={ 73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,72,71,69,68,66,64,62,58,54,49,41,33,25,17,13,15,21,30,39,47,54,59,62,64,66,67,68,69,70,71,71,72,72,72,71,71,70,69,68,67,66,65,63,62,60,56,51,45,37,29,22,18,18,22,28,33,35,36,35,32,26,20,15,12,15,20,26,31,35,37,36,34,30,25,22,22,27,33,41,48,55,60,63,66,67,68,69,70,71,72,72,73,73,73,73,73,73,72,71,70,69,67,65,63,60,55,49,40,30,21,13, 7,10,17,27,36,45,52,56,59,60,61,62,62,62,62,61,61,59,57,53,47,40,32,24,18,15,18,23,28,32,33,31,28,23,16,10, 6, 8,13,20,27,31,32,31,28,22,15,10, 6,10,16,23,30,34,36,36,34,29,24,20,19,24,30,37,44,51,56,59,61,62,63,64,64,64,65,64,64,62,60,57,53,48,43,38,36,39,43,49,54,59,63,66,68,69,70,71,72,72,73,73,73,73,73,73,73,73,73,73,73,73,73,73 };
float f1[]={ 55,58,60,62,64,66,67,68,68,69,69,70,71,72,72,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,73,72,72,72,72,72,72,72,72,72,72,72,72,72,72,72,73,73,73,72,72,72,72,71,71,71,71,71,71,70,70,69,68,67,66,64,63,61,60,59,57,55,52,49,46,43,40,37,35,34,33,32,32,33,34,36,37,39,41,43,45,47,50,52,55,57,59,60,61,61,62,62,62,62,61,61,60,58,57,55,53,51,49,48,46,44,42,40,38,35,32,30,29,28,27,27,26,26,26,25,25,24,23,23,23,24,24,25,25,26,26,26,27,28,29,31,33,35,38,40,41,43,44,46,48,50,53,55,57,59,60,61,62,63,64,64,65,65,64,63,61,59,57,54,52,50,47,45,42,39,37,34,32,31,30,30,30,31,32,34,36,37,39,40,41,42,43,44,44,44,44,43,42,41,40,38,36,34,32,30,28,26,25,24,23,22,21,20,18,17,17,17,17,18,18,18,18,18,18,18,18,18,18,19,19,19,19,20,20,21,23,24,25,26,26,26,27,28,29,31,34,36,37,38,40,41,43,45,47,48,49,50,51,51,51,50,49,49,48,48,47,47,47,47,47,47,48,60 };
const int n=sizeof(f0)/sizeof(f0[0]);

所有值都需要转换：

f0[i] = 73.0-f0[i];
f1[i] = 73.0-f1[i];

要从图像偏移...f0 是原始红色信号，f1 是失真的黄色信号。

这是我使用一阶 FIR 滤波器得到的最接近的结果：

上半部分是使用的 FIR 滤波器权重图（可通过鼠标编辑，因此 FIR 是手绘的，以便快速找到权重...）。下面是信号图：

• 红色原始（理想）信号f0

• 深绿色测量信号f1

• 浅绿色 FIR 滤波理想信号f2

  FIR 滤波器只是卷积，其中零偏移元素是最后一个，这里是权重值：

float fir[35] = { 0.0007506932, 0.0007506932, 0.0007506932, 0.0007506932, 0.0007506932, 0.0007506932, 0.0007506932, 0.0007506932, 0.0007506932, 0.003784107, 0.007575874, 0.01060929, 0.01591776, 0.02198459, 0.03032647, 0.04170178, 0.05686884, 0.06445237, 0.06900249, 0.07203591, 0.07203591, 0.0705192, 0.06900249, 0.06672744, 0.06217732, 0.05611049, 0.04928531, 0.04170178, 0.03335989, 0.02653471, 0.02046788, 0.01515941, 0.009850933, 0.005300813, 0.0007506932 };

所以要么有更高程度的 FIR，要么需要对权重进行更多调整。无论如何，这对于去卷积和/或拟合应该足够了......顺便说一句，FIR 滤波器的完成如下：

const int fir_n=35;     // size (exposure time) [samples]
const int fir_x=fir_n-1;    // zero offset element [samples]
int x,y,i,j,ii;
float a,f2[n];
for (i=0;i<n;i++)
 for (f2[i]=0.0,ii=i-fir_x,j=0;j<fir_n;j++,ii++)
  if ((ii>=0)&&(ii<n)) f2[i]+=fir[j]*f0[ii];

Answer 1

在我看来，您应该从一个非常简单的系统识别和连续的信号重建开始。另外，我建议首先在 Matlab（商业许可）或 Octave（免费→https://www.gnu.org/software/octave/download.html）等数学原型工具中实现您的算法。无论您使用什么库，这些工具都提供了 Pascal 或 Java 等编程语言无法提供的轻松信号处理。用 Matlab 或 Octave 成功设计算法后，再考虑如何用 Pascal 实现它。

让我们假设电子管的行为可以通过线性时不变系统（例如线性低通滤波器）来表征。这绝不是保证，而是一种值得的方法（至少在它失败之前:)）。对非线性和/或时变系统采用相同的方法会变得相当复杂，我想你需要专业的帮助才能做到这一点。

如果我正确理解了您的描述，您可以访问电子管的输入和输出信号。如果我错了，您不知道输入信号，您可以先应用一些您知道其特性的校准信号，然后记录输出信号。了解输入和输出信号是以下方法的先决条件。如果没有这两个信号，您将无法估计电子管的脉冲响应 h。在计算出 h 的近似值后，我们可以设计一个称为 ge 的逆滤波器，并最终从输出信号重构输入。

这是输入信号 x[n] 通过您的管子的信号流，由 h 描述，产生输出信号 y[n]。取 y[n] 并应用由 ge 描述的逆滤波操作，我们得到 xr[n]

x[n] →| h | → y[n] → |通用电气 | → xr[n]

取一个长度为 N 的输入向量 x 和一个对应的长度相同的 y 向量。 现在，您将输出 y 表示为输入卷积矩阵 X（参见下面的代码以了解其实现）与系统的未知脉冲响应的卷积，即

y = X * h

向量和矩阵大小为 y = N x 1、X = N x N 和 h = N x 1 您可以通过计算

来计算脉冲响应 he 的最小二乘近似

他 = inv(X'*X)*X' * y

其中 X' 描述了转置，inv() 描述了 X 的逆矩阵。他表示我们通过一维反卷积获得的已识别电子管脉冲响应的列向量。您可以通过计算估计系统的输出来检查识别的效果，

ye = X * he

并通过比较 ye 和 y。现在，我们尝试从 y 和 he 重构 x。重构后的输入向量 xr 由

计算

xr = Ge * y

其中 Ge = inv(He) 并且 He 是 he 的 N x N 卷积矩阵。 这是一些八度代码。将这两个函数复制到它们自己的专用文件（reconstruct.m 和 getConvolutionMatrix.m）中，然后在 Octave 命令行中键入“reconstruct.m”以检查示例的输出。请注意，示例代码仅适用于奇数长度的向量（N 为奇数）。玩弄你的向量的大小 N。这可能有助于近似精度。

function [Ge] = reconstruct ()
    x = [1 2 3]';            # input signal
    # h = [1 3 2]';          # unknown impulse response
    y = [5 11 13]';          # output signal
    y = y + 0.001*randn(length(y),1)  # add noise to output signal 

    Xm = getConvolutionMatrix(x)
    Xps = inv(Xm'*Xm)*Xm';
    he = Xps * y
    He = getConvolutionMatrix(he);
    Ge = inv(He);
    # reconstructed impulse signal
    xr = Ge*y
endfunction

function [mH] = getConvolutionMatrix(h)
    h = h(:)';
    hlen = length(h);
    Nc = (hlen-1)/2;

    mH= zeros(hlen, hlen);
    hp = [zeros(1,Nc) h zeros(1,Nc)];

    for c=1:hlen
      for r=1:hlen
        mH(r,c) = hp(r+hlen-c);
      end
    end
endfunction