信号增强算法 [英] Signal enhancing algorithm

问题描述

我需要一种算法(最好是一种类似于Pascal的语言，但是最终并不重要)，它将使信号"成为可能. (实际上是一系列数据点)在左侧看起来像在右侧.

I need an algorithm (preferable in a Pascal-like language, but it the end doesn't really matter) that will make the "signal" (actually a series of data points) in left look like the one in right.

信号来源:
该信号是由机器产生的.为了简化说明，该机器正在测量流过透明管的液体的密度.因此，该信号类似于电信号(音频/射频). 数据点可能看起来像这样:[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]

Signal origin:
The signal is generated by a machine. Oversimplifying the explanation, the machine is measuring the density of a liquid flowing through a transparent tube. So, the signal is nothing similar to an electrical signal (audio/radio frequency). The data points could look like this: [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]

我需要什么:
我需要准确检测峰值".为此，我已经有一段代码，但是它不能处理不良信号，如下图所示.
我认为我们可以看到这一点，因为它是偶然通过低通滤波器的信号，现在我想恢复它.

What I need:
I need to accurately detect the 'peaks'. For this, I already have a piece of code but it is not working on poor signals as the one shown in the image below.
I think we can see this, as a signal that was accidentally passed through a low-pass filter, and now I want to restore it.

注意:
有4个信号，但是它们是分开的，因此可以分别进行分析.因此，仅考虑如何处理其中之一就足够了.

Notes:
There are 4 signals, but they are separated so they can be analyzed individually. So, it is enough to think about how to process just one of them.

在一个峰值之后，如果信号下降得不够快，我们可以认为存在多个峰值(您最好在系列末尾的红色"信号中看到该峰值).

After a peak, if the signal is not coming down fast enough we can consider that there are multiple peaks (you can best see that in the 'red' signal at the end of the series).

优点是整个系列都可用(因此信号不是实时的，它已经存储在文件中了！)

The advantage is that the whole series is available (so the signal is not in real time, it is already stored on file)!

[由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是失真的黄色信号.

To offset back from image... f0 is the original Red signal and f1 is the distorted Yellow one.

这是我与一阶FIR滤波器最接近的一面:

This is the closest I get to with first order FIR filter:

上半部分是使用的FIR滤波器权重的图(可通过鼠标进行编辑，因此可以手工绘制FIR以快速找到权重...).下面是信号图:

The upper half is plot of used FIR filter weights (editable by mouse so the FIR is hand drawed for fast weights find...). Bellow are the signal plots:

• 红色原始(理想)信号f0

深绿色测量信号f1

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

Light Green FIR filtered Ideal signal f2

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

The FIR filter is just convolution where zero offset element is the last and here the weight values:

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滤波器按如下方式完成:

So either there is also some higher degree of FIR or the weights need to be tweaked a bit more. Anyway this should be enough for the deconvolution and or fitting ... btw the FIR filter is done as follows:

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];

推荐答案

我认为，您应该从一个非常简单的系统标识和连续的信号重建开始.另外，我建议您先在数学原型工具(例如Matlab(商业许可)或Octave(免费→ https://www.gnu.org/software/octave/download.html ).这些工具提供了轻松的信号处理能力，无论您使用哪种库，都无法像Pascal或Java那样提供编程语言.使用Matlab或Octave成功设计算法后，请考虑如何使用Pascal实施该算法.

In my opinion, you should start with a pretty simple system identification and consecutive signal reconstruction. Also, I would recommend to implement your algorithm first in a mathematical prototyping tool like Matlab (commerical license) or Octave (free → https://www.gnu.org/software/octave/download.html). These tools provide an ease of signal processing no programming language like Pascal or Java could ever offer, no matter what library you use. After you successfully designed your algorithm with Matlab or Octave, then think about how to implement it with Pascal.

让我们假设电子管的行为可以通过线性时不变系统(例如线性低通滤波器)来表征.这绝不是可以保证的，而是值得采用的方法(至少直到失败为止:)).对于非线性和/或时变系统，采用相同的方法非常困难，我想您将需要专业的帮助.

Lets assume the behaviour of the tube can be characterized by a linear time-invariant system (e.g. a linear lowpass filter). This is by no means guaranteed but a worthwhile approach (at least until it fails:) ). Following the same approach for non-linear and/or time-variant systems becomes pretty involved and you will need professional help to do this I would imagine.

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

If I understood your description correctly, you have access to both the input and output signals of the tube. If I am wrong and you don’t know the input signal, you might be able to apply some calibration signal first, whose characteristics you know, and record the output signal. Knowing input and output signals is a prerequisite for the following approach. Without both signals you can not approximate the impulse response h of the tube. After calculating an approximation of h, we can design an inverse filer called ge and eventually reconstruct the input from the output signal.

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

Here is the signal flow of a input signal x[n] passing your tube described by h producing the output signal y[n]. Taking y[n] and applying an inverse filtering operation described by ge we obtain xr[n]

x [n]→| h | →y [n]→| ge | →xr [n]

取长度为N的输入向量x和长度相同的y的对应向量. 现在，您将输出y表示为输入卷积矩阵X的卷积(有关实现，请参见下面的代码)，其中包含系统未知的脉冲响应，即

Take an input vector x of length N and a corresponding vector of y of the same length. Now you express the output y as a convolution of an input convolution matrix X (see the code below for its implementation) with the unknown impulse response of your system, i.e.

y = X * h

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

with the vector and matrix sizes y = N x 1, X = N x N and h = N x 1 You can calculate a least-squares approximation of the impulse response he, by calculating

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

其中X'描述了X的转置和inv()的矩阵逆.他表示已确定的您的管的冲激响应的列向量，该列向量是通过 1-D反卷积获得的.您可以通过计算估算系统的输出来检查标识的工作情况，

where X' describes the transpose and inv() the matrix inverse of X. he represents a column vector of the identified impulse response of your tube which we obtained via a 1-D deconvolution. You can check how well the identification worked by calculating the output of your estimated system,

ye = X *他

并通过比较ye和y.现在，我们尝试从y和他重建x.重建的输入向量xr由

and by comparing ye and y. Now, we try to reconstruct x from y and he. The reconstructed input vector xr is calculated by

xr = Ge * y

其中Ge = inv(He)，He是He的N x N卷积矩阵. 这是一些倍频程代码.将这两个函数复制到各自的专用文件中(reconstruct.m和getConvolutionMatrix.m)，然后在Octave命令行中键入"reconstruct.m"以检查示例的输出.请注意，示例代码仅适用于奇数长度的向量(N为奇数).播放N个向量的大小.这可能有助于近似精度.

where Ge = inv(He) and He is the N x N convolution matrix of he. Here is some Octave code. Copy both functions in their own dedicated file (reconstruct.m and getConvolutionMatrix.m) and type "reconstruct.m" into the Octave command line to examine the outputs of the example. Please note the sample code works only for vector of odd length (N is odd). Play around with the size N of your vectors. This might help the approximation accuracy.

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

这篇关于信号增强算法的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！