Output of variable known, impulse response is known, how to get input signal in matlab?

Hello all, I have an oxygen sensor that has reported back oxygen
concentrations while moving in a yo-yo pattern through the ocean. 

I have a vector of oxygen concentrations (about 10000 entries or more) ,
and a vector of corresponding time intervals (approximately constant ~
every 2 seconds). I consider my o2 concentration the output signal y(t).

I know the sensor responds exponentially, so we are assuming an impulse
response h(t)=[1-(e^-t/tau)] where tau is our time constant ~ 24seconds.

Can anyone out there suggest how I would go about solving for the initial
input signal x(t)? I assume this can be done in matlab, perhaps the signal
processing toolbox.

Any help, references, comments are welcome.

cheers,

- Charlie

Charlie wrote:
> Hello all, I have an oxygen sensor that has reported back oxygen
> concentrations while moving in a yo-yo pattern through the ocean.
>
> I have a vector of oxygen concentrations (about 10000 entries or more) ,
> and a vector of corresponding time intervals (approximately constant ~
> every 2 seconds). I consider my o2 concentration the output signal y(t).
>
> I know the sensor responds exponentially, so we are assuming an impulse
> response h(t)=[1-(e^-t/tau)] where tau is our time constant ~ 24seconds.
>
> Can anyone out there suggest how I would go about solving for the initial
> input signal x(t)? I assume this can be done in matlab, perhaps the signal
> processing toolbox.
>
> Any help, references, comments are welcome.

Hello Charlie,

welcome back. You increased your sampling rate from 4s to 2s, that is
good (could be better, though). Depending on your yo-yo rate, you
might still have pretty severe aliasing problems that could render
your measurements useless.

Your first task is to re-sample the measurements to a constant
sampling period. Then try the inverse filter that I mentioned in the
last thread and see if the output is close to what you expect.

Regards,
Andor

"charliebishop" <charliebishop@gmail.com> wrote in message 
news:2O-dnbMylYomm3_YnZ2dnUVZ_vCknZ2d@giganews.com...
> Hello all, I have an oxygen sensor that has reported back oxygen
> concentrations while moving in a yo-yo pattern through the ocean.
>
> I have a vector of oxygen concentrations (about 10000 entries or more) ,
> and a vector of corresponding time intervals (approximately constant ~
> every 2 seconds). I consider my o2 concentration the output signal y(t).
>
> I know the sensor responds exponentially, so we are assuming an impulse
> response h(t)=[1-(e^-t/tau)] where tau is our time constant ~ 24seconds.
>
> Can anyone out there suggest how I would go about solving for the initial
> input signal x(t)? I assume this can be done in matlab, perhaps the signal
> processing toolbox.
>
> Any help, references, comments are welcome.
>
> cheers,
>
> - Charlie

Charlie,

You still haven't given us critical information that might help in answering 
your question.  Here would be a few:

- what is the rate of descent of the sensor?

- what is the rate of ascent of the sensor?

- what is the maximum expected rate of change of oxygen concentration as a 
function of depth?

- given the above, what is the maximum expected rate of change of oxygen 
concentration per unit time?  (just as a check on the numbers).

Then, I would be tempted to construct a system model with the motion (depth) 
as a forcing function and then solve for an input that creates the known 
output - a prediction system of sorts.  I think I already sent you a diagram 
of something like that.

The fact remains, if information has been destroyed by the sampling rate, 
there's no getting it back!
The fact remains that since the system is a simple integrator or averager or 
lowpass filter then the inverse system is a differentiator.  Differentiators 
are "noisy".  That is, their output will be noisy looking if there are any 
jumps in their input.

Usually when doing a differentiator as processor I have had to do this:

1) Lowpass the data as much as possible ahead of the differentiator.
2) Lowpass the data as much as possible after the differentiator.
In linear system theory it should be possible to interchange the processing 
boxes and do all the filtering before or after.  However, I find it nice to 
be able to *observe* the input to the differentiatorinput filter (the raw 
data), the input to the differentiator, the output of the differentiator and 
the filtered output.  That way each of the two filters can be adjusted to 
suit your needs better I think.  That is, it's easier to observe the input 
to the differentiator and decide when the input has been filtered "too 
much".  Ditto at the output.
Then I use ramp inputs or triangular wave inputs to get the transient 
response I want.  The output from a triangular wave should be a square wave 
with exponentially rising and falling edges.  The tradeoff is:
- do the square wave transitions reach steady state soon enough?
- is the noise level small enough?
The former is easy to do with synthetic data.
The latter probably needs real data to adjust.

Fred 

Okay, before answering these questions, allow me to give a little
better explanation of what it is I'm trying to do.

We have an autonomous underwater vehcile which profiles the ocean,
surface to 200m depth in a yo-yo pattern.. each transect is about 8km
long (several hundred yo's in some cases)

The instrument is measuring conductivitiy, temp (From ctd) , pressure,
vehicle speed, oxygen concentration, etc.

The problem is because there is a strong vertical gradient in oxygen
in our study area, when we plot the downcast of the vehicle with the
upcast of the vehicle, the profiles for oxygen do not match up (But
they should) ... the profiles do match for our CTD, just not our
oxygen sensor because of its slower sampling rate.

So I dont necessarily want to find out the original input signal, I
just want to adjust the data so that the upcast data matches the
downcast.. because right now we don't know which is right... the
reason I want to compare them is because I can see an interesting
phenomena in the downcast data, that does not appear in the
corresponding upcast... but then again, because of the slow response
time, when I look at the interesting data I'm not even sure what depth
it corresponds to!

so basically I want to adjust my oxygen concentrations so that they
correspond to the depth the vehicle is located at, at any one time.

The descent and ascent of my vehicle are appoximately the same, about
20cm/s

the max oxygen concentration (at the surface) is about 15mg/L , its
minimum is about 8mg/L.

Any thoughts, comments, criticisms?

- Charlie



On Feb 25, 8:12 pm, "Fred Marshall" <fmarshallx@remove_the_x.acm.org>
wrote:
> "charliebishop" <charliebis...@gmail.com> wrote in message
>
> news:2O-dnbMylYomm3_YnZ2dnUVZ_vCknZ2d@giganews.com...
>
>
>
>
>
>
>
> > Hello all, I have an oxygen sensor that has reported back oxygen
> > concentrations while moving in a yo-yo pattern through the ocean.
>
> > I have a vector of oxygen concentrations (about 10000 entries or more) ,
> > and a vector of corresponding time intervals (approximately constant ~
> > every 2 seconds). I consider my o2 concentration the output signal y(t).
>
> > I know the sensor responds exponentially, so we are assuming an impulse
> > response h(t)=[1-(e^-t/tau)] where tau is our time constant ~ 24seconds.
>
> > Can anyone out there suggest how I would go about solving for the initial
> > input signal x(t)? I assume this can be done in matlab, perhaps the signal
> > processing toolbox.
>
> > Any help, references, comments are welcome.
>
> > cheers,
>
> > - Charlie
>
> Charlie,
>
> You still haven't given us critical information that might help in answering
> your question.  Here would be a few:
>
> - what is the rate of descent of the sensor?
>
> - what is the rate of ascent of the sensor?
>
> - what is the maximum expected rate of change of oxygen concentration as a
> function of depth?
>
> - given the above, what is the maximum expected rate of change of oxygen
> concentration per unit time?  (just as a check on the numbers).
>
> Then, I would be tempted to construct a system model with the motion (depth)
> as a forcing function and then solve for an input that creates the known
> output - a prediction system of sorts.  I think I already sent you a diagram
> of something like that.
>
> The fact remains, if information has been destroyed by the sampling rate,
> there's no getting it back!
> The fact remains that since the system is a simple integrator or averager or
> lowpass filter then the inverse system is a differentiator.  Differentiators
> are "noisy".  That is, their output will be noisy looking if there are any
> jumps in their input.
>
> Usually when doing a differentiator as processor I have had to do this:
>
> 1) Lowpass the data as much as possible ahead of the differentiator.
> 2) Lowpass the data as much as possible after the differentiator.
> In linear system theory it should be possible to interchange the processing
> boxes and do all the filtering before or after.  However, I find it nice to
> be able to *observe* the input to the differentiatorinput filter (the raw
> data), the input to the differentiator, the output of the differentiator and
> the filtered output.  That way each of the two filters can be adjusted to
> suit your needs better I think.  That is, it's easier to observe the input
> to the differentiator and decide when the input has been filtered "too
> much".  Ditto at the output.
> Then I use ramp inputs or triangular wave inputs to get the transient
> response I want.  The output from a triangular wave should be a square wave
> with exponentially rising and falling edges.  The tradeoff is:
> - do the square wave transitions reach steady state soon enough?
> - is the noise level small enough?
> The former is easy to do with synthetic data.
> The latter probably needs real data to adjust.
>
> Fred

Also, please see the following image for an example of the difference
between the upcast and downcast data for one particular transect.
Upcast and downcast data should be the same because the water column
being sampled by each is more or less the same.

http://www.physics.mun.ca/~cbishop/images/anomaly1.jpg

cheers,

- Charlie

Also, please see the following image for an example of the difference
between the upcast and downcast data for one particular transect.
Upcast and downcast data should be the same because the water column
being sampled by each is more or less the same.

http://www.physics.mun.ca/~cbishop/images/anomaly1.jpg

cheers,

- Charlie

Charlie wrote:
> Okay, before answering these questions, allow me to give a little
> better explanation of what it is I'm trying to do.
> 
> We have an autonomous underwater vehcile which profiles the ocean,
> surface to 200m depth in a yo-yo pattern.. each transect is about 8km
> long (several hundred yo's in some cases)
> 
> The instrument is measuring conductivitiy, temp (From ctd) , pressure,
> vehicle speed, oxygen concentration, etc.
> 
> The problem is because there is a strong vertical gradient in oxygen
> in our study area, when we plot the downcast of the vehicle with the
> upcast of the vehicle, the profiles for oxygen do not match up (But
> they should) ... the profiles do match for our CTD, just not our
> oxygen sensor because of its slower sampling rate.
> 
> So I dont necessarily want to find out the original input signal, I
> just want to adjust the data so that the upcast data matches the
> downcast.. because right now we don't know which is right... the
> reason I want to compare them is because I can see an interesting
> phenomena in the downcast data, that does not appear in the
> corresponding upcast... but then again, because of the slow response
> time, when I look at the interesting data I'm not even sure what depth
> it corresponds to!

The lag makes both upward and downward readings wrong. Aside from 
correcting for the DO sensor's lag (which you have learned by now how to 
do) the best you can do is smooth the data (which adds a known lag you 
can correct for) then average readings for the same depth from upward 
and downward passes. The greatest error will be at mid depth, where the 
two readings will have been taken farthest apart.

> so basically I want to adjust my oxygen concentrations so that they
> correspond to the depth the vehicle is located at, at any one time.
> 
> The descent and ascent of my vehicle are appoximately the same, about
> 20cm/s
> 
> the max oxygen concentration (at the surface) is about 15mg/L , its
> minimum is about 8mg/L.
> 
> Any thoughts, comments, criticisms?

Why the resistance to removing the lag by deconvolution?

Jerry
-- 
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Charlie wrote:
> Also, please see the following image for an example of the difference
> between the upcast and downcast data for one particular transect.
> Upcast and downcast data should be the same because the water column
> being sampled by each is more or less the same.

No upcast measurements should represent conditions at a greater depth 
and downcast measurements should represent measurements at a lesser 
depth, which is what I seem to see. Other parameters exhibit hysteresis too.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

On 7 Mar, 01:11, "Charlie" <charliebis...@gmail.com> wrote:
> Also, please see the following image for an example of the difference
> between the upcast and downcast data for one particular transect.
> Upcast and downcast data should be the same because the water column
> being sampled by each is more or less the same.
>
> http://www.physics.mun.ca/~cbishop/images/anomaly1.jpg

Interesting profiles indeed.

The pattern is too consistent between the two deep
scans to be a coincidence. The time axes are impossible
to make sense of. How long time did these deep scans
take? How many samples in the down and upscans?
What is the horizontal speed of the towfish?

As for how to continue, it depends on exactly how the
oxigen sensor works. Is the function you mention in
your first post the impulse response or the step
response? From plotting it, I assume it is the step
response, but I don't know these sorts of systems
so I may very well be wrong.

Assuming this really is the step response, find the
*impulse* response as the time derivative of the step
response. It should come out something like

h(t) = d/dt[1-exp(-t/tau)]=1/tau*exp(-t/tau)

but don't take my word for it! Arithmetics off the
top of my head has never been my force; a night
of insomnia certainly doesn't help in that
department.

Next, sample h(t) and compute the z transform, ZT:

h[n] = 1/tau*exp(-nT/tau)= [1/tau*exp(-T/tau)]^n

where T is the sampling period. The above is a
standard form where the ZT is known to be

H(z) = 1/(1-[1/tau*exp(-T/tau)]*z^-1)

This is a 1st order IIR filter, which has a FIR
inverse h'[n] as

h'[n] =[1, -1/tau*exp(-T/tau)]

Plot your oxygen data as a time series, convolve
with h'[n] and see if there are improvements.

Now, the above comes with all the inherent
guarantees of USENT (i.e. none whatsoever), my
insomnia, aritmetics and guesswork spicing up the
mix. Do note that I haven't checked, only assumed,
that the Nyquist criterion holds and that the pole
of h[n] actually is inside the unit circle.

If there still are discrepancies between the downcast
and upcast after filtering, keep in mind that the
downcast is likely to go through undisturbed water,
while the upcast is likely to go through turbulent
water. The cable between the tow-fish and the vessel
has at least disturbed the water; maybe even the
downcast has set up some turbulence that affects the
upcast. You might want to have a chat with the
marine engineers and oceanographers to find out
the likely effects of such possible cable-induced
turbulence.

Rune

Jerry Avins wrote:
> Charlie wrote:
>> Also, please see the following image for an example of the difference
>> between the upcast and downcast data for one particular transect.
>> Upcast and downcast data should be the same because the water column
>> being sampled by each is more or less the same.
> 
   No. Upcast measurements should represent conditions at a greater depth
> and downcast measurements should represent measurements at a lesser 
> depth, which is what I seem to see. Other parameters exhibit hysteresis 
> too.
> 
> Jerry


-- 
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯