time series - Signal enhancing algorithm -

i need algorithm (preferable in pascal-like language, end doesn't matter) make "signal" (actually series of data points) in left 1 in right.

signal origin:
signal generated machine. oversimplifying explanation, machine measuring density of liquid flowing through transparent tube. so, signal nothing similar electrical signal (audio/radio frequency). data points 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 need:
need accurately detect 'peaks'. this, have piece of code not working on poor signals 1 shown in image below.
think can see this, signal accidentally passed through low-pass filter, , want restore it.

notes:
there 4 signals, separated can analyzed individually. so, enough think how process 1 of them.

after peak, if signal not coming down fast enough can consider there multiple peaks (you can best see in 'red' signal @ end of series).

the advantage whole series available (so signal not in real time, stored on file)!

[edit spektre] extracted red sample points image

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

all values need transformed:

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

to offset image... f0 original red signal , f1 distorted yellow one.

this closest first order fir filter:

the upper half plot of used fir filter weights (editable mouse fir hand drawed fast weights find...). bellow signal plots:

• red original (ideal) signal f0
• dark green measured signal f1

• light green fir filtered ideal signal f2

  the fir filter convolution 0 offset element last , here 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 }; 

so either there higher degree of fir or weights need tweaked bit more. anyway should enough deconvolution , or fitting ... btw fir filter done follows:

const int fir_n=35;     // size (exposure time) [samples] const int fir_x=fir_n-1;    // 0 offset element [samples] int x,y,i,j,ii; float a,f2[n]; (i=0;i<n;i++)  (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]; 

in opinion, should start pretty simple system identification , consecutive signal reconstruction. also, recommend implement algorithm first in mathematical prototyping tool matlab (commerical license) or octave (free → https://www.gnu.org/software/octave/download.html). these tools provide ease of signal processing no programming language pascal or java ever offer, no matter library use. after designed algorithm matlab or octave, think how implement pascal.

lets assume behaviour of tube can characterized linear time-invariant system (e.g. linear lowpass filter). no means guaranteed worthwhile approach (at least until fails:) ). following same approach non-linear and/or time-variant systems becomes pretty involved , need professional imagine.

if understood description correctly, have access both input , output signals of tube. if wrong , don’t know input signal, might able apply calibration signal first, characteristics know, , record output signal. knowing input , output signals prerequisite following approach. without both signals can not approximate impulse response h of tube. after calculating approximation of h, can design inverse filer called ge , reconstruct input output signal.

here signal flow of input signal x[n] passing tube described h producing output signal y[n]. taking y[n] , applying inverse filtering operation described ge obtain xr[n]

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

take input vector x of length n , corresponding vector of y of same length. express output y convolution of input convolution matrix x (see code below implementation) unknown impulse response of system, i.e.

y = x * h

with vector , matrix sizes y = n x 1, x = n x n , h = n x 1 can calculate least-squares approximation of impulse response he, calculating

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

where x' describes transpose , inv() matrix inverse of x. represents column vector of identified impulse response of tube obtained via 1-d deconvolution. can check how identification worked calculating output of estimated system,

ye = x * he

and comparing ye , y. now, try reconstruct x y , he. reconstructed input vector xr calculated

xr = ge * y

where ge = inv(he) , n x n convolution matrix of he. here octave code. copy both functions in own dedicated file (reconstruct.m , getconvolutionmatrix.m) , type "reconstruct.m" octave command line examine outputs of example. please note sample code works vector of odd length (n odd). play around size n of vectors. might 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 output signal       xm = getconvolutionmatrix(x)     xps = inv(xm'*xm)*xm';     = xps * y     = 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)];      c=1:hlen       r=1:hlen         mh(r,c) = hp(r+hlen-c);       end     end endfunction