Am 1/28/2013 12:15 PM, schrieb Tim Wescott:
> On Mon, 28 Jan 2013 11:46:02 -0800, Peter Mairhofer wrote:
> 
> [...]
>> My observation from simple MATLAB experiments is:
>>
>> 1) I can use any subset of the equations and get a solution. This is in
>> one sense obvious because it is LS. On the other hand, it contradicts my
>> intuition. Suppose my Nyquist rate for input x(t) and the unknown system
>> h(t) is 10kHz and h(t) operates (amongst others) also at 4 kHz. If I
>> sample the output y(t) at 5kHz, how should the reconstruction algorithm
>> know about the 4kHz range?
> 
> Oy.  Nyquist didn't say that.  In fact, your above statement makes almost 
> no sense.  http://www.wescottdesign.com/articles/Sampling/sampling.pdf.
> 
> You state elsewhere that you are modeling the system h(t) as something 
> that can be reasonably represented as a FIR filter.  In that case, the 
> assumption that h(t) can be reasonably represented as a FIR filter is 
> only valid if you do your measurements at the same sampling rate as your 
> FIR filter is going to operate.

Maybe I should state the setup more precisly:
- I assume an ideal (mathematical) setup (exact pointwise samples, ideal
low-pass filters, infinite long signals)
- I assume no noise for now
- I assume perfect linearity
- I assume that my input signal x(t) perfectly bandlimited white noise
between 0-5kHz and is zero >5kHz.

The filter h(t) will operate on x(t), modifying any frequency content
between 0 and 5kHz (it can't do anything more because it is a perfect
linear system).

Because of the assumptions described above, I can convert the whole
setup to a discrete-time setup: x[n] represents x(t), y[n] y(t) where
the signals are sampled at n*Ts (Ts=1/10kHz). h(t) can be converted to
an IIR filter via bilinear transform (operating on x[n] the same way as
on x(t) between 1-5kHz).

The only non-ideal assumption: The IIR equivalent of h(t) is
approximated as an FIR filter h[n] with Q taps.

Assuming that h[n] models h(t) good enough, N samples of x[n] and y[n]
should recover the coefficients of h[n]:

h = X^+ * y

where (.)^+ denotes the Moore-Penrose pseudo inverse, X is a matrix
containg the tap-inputs of x[n] and y[n] contains N samples.

Stated differently: When I feed x[n] to my IIR equivalent of h(t) and
h[n] via y=X*h, both should produce a similar y[n] (as similar as h[n]
approximates the IIR equivalent of h(t)).

Now to what I meant in my original question: Is it possible to use just
any second or forth sample in y[n] and still recover h[n] perfectly?
Note that "every second sample" corresponds to a subsampling of 2 in the
ideal setup.

The answer is yes (in theory) because it is LS (which is solveable as
long as there are Q lin. independent rows).

What I did next is to sum over M consecutive samples of y[n] and take
the sum. This "sum-and-dump" corresponds to subsampling of factor M with
pre-lowpass filtering (in the ideal setup). In order to hold the
equality in y=X*h, I need to apply the summation in X as well. This
approach still recovered the h[n] under perfect numerical conditions
although I destroyed much frequency content because of the summation.
However, because of the summation in x[n] over M samples, the columns
become more and more orthogonal and the system therefore numerically
instable.

Please regard this as a tough-experiment rather than a question
regarding practical SysId.

Peter

On Mon, 28 Jan 2013 11:46:02 -0800, Peter Mairhofer wrote:

> Hi,
> 
> Thank you (and all) for the discussion. I searched some typical books in
> adaptive filter theory but this questions is never addressed.
> 
> Am 1/28/2013 9:20 AM, schrieb Tim Wescott:
>> On Mon, 28 Jan 2013 00:09:51 -0800, Peter Mairhofer wrote:
>> 
>>> Hi,
>>>
>>> Is there a theoretical minimum for the sampling rate of the output
>>> y(t) in a SysId setup?
>> 
>> The dangerous word in your question is "theoretical".  If you never
>> need a usable answer, then see Vladimir's response, because yes, in
>> theory you could do this at just about any sampling rate.
> 
> Ok, thank you. Acutally, I am interested in both cases - theoretical and
> practical to understand my observations.
> 
> For the theoretical side, think about that: Is it a fundamental barrier
> or could it beaten by a method which is not yet discovered? E.g., we
> know now that the Nyquist rate is not the fundamental barrier for the
> sampling rate, it's the information content (c.f. Compressive Sampling).
> It is just a matter if it is possible to develop an appropriate sampling
> scheme exploiting a specific structure.
> 
> So right now it sounds to me that there is no such fundamental limit
> which has been proven. So it could be possible to discover completely
> new ways of SysId, sampling the output at a lower rate which could even
> be physically reliable.

You need to be very very careful about making any claims about physical 
reliability.

> My observation from simple MATLAB experiments is:
> 
> 1) I can use any subset of the equations and get a solution. This is in
> one sense obvious because it is LS. On the other hand, it contradicts my
> intuition. Suppose my Nyquist rate for input x(t) and the unknown system
> h(t) is 10kHz and h(t) operates (amongst others) also at 4 kHz. If I
> sample the output y(t) at 5kHz, how should the reconstruction algorithm
> know about the 4kHz range?

Oy.  Nyquist didn't say that.  In fact, your above statement makes almost 
no sense.  http://www.wescottdesign.com/articles/Sampling/sampling.pdf.

You state elsewhere that you are modeling the system h(t) as something 
that can be reasonably represented as a FIR filter.  In that case, the 
assumption that h(t) can be reasonably represented as a FIR filter is 
only valid if you do your measurements at the same sampling rate as your 
FIR filter is going to operate.

You're basically sitting on a steam roller here, driving over too many 
mathematical subtleties to list.

> 2) Using a lowpass filter to supress aliasing before sampling still
> gives a valid result when the numbers are absolutely accurate. Minor
> perturbations however, give absolutely wrong results.

Without knowing how you are applying the lowpass filter, it's hard to 
comment on the validity of that assertion.  In general, if you excite the 
system with a signal x(t) and look at the output y(t), and get the result 
h(t), then if you turn around and filter both x(t) and y(t) by the same 
identical filter to get x'(t) and y'(t), then using those two filtered 
signals in your system identification should cough up the same h(t).

But, you will be more sensitive to noise at high frequencies.

> My conclusion therefore is: The (implicit) aliasing when just using
> *perfect* pointwise values of y(t) does the job. Low-pass filtering
> removes the vital information from the system, rendering the
> reconstruction very unstable.

I think you need to stop trying to think of this in terms of aliasing.  
The frequency domain is a great place to work, but it is just a 
mathematical shortcut to solving certain problems in linear systems 
design.  This problem isn't one of those, so trying to think in frequency 
domain terms will just confuse you.

It's a time domain problem.  Keep it there.

> (2) even gives correct results under perfect numerical conditions. My
> mathematical interpretations for this is: The columns in the matrix
> become nearly linearly dependent because of the "smoothing".

This could well be, depending on how much smoothing you're doing.

>>> I know that the input signal x(t) must contain all frequency
>>> components ("persistently exciting") of the unknown filter.
>> 
>> You know wrong. The input signal must contain _enough_ frequency
>> components to persistently excite the unknown filter. [...]
> 
> Sure. Sorry, that's what I meant.
> 
>>> When using a simple LS method to solve for the coefficients, the
>>> overdetermined system of equations suggests that it is possible to
>>> drop rows (as long as there are enough linearly independent) and thus
>>> (implicitely) subsample y(t).
>> 
>> Yes, that is the suggestion.  And if you're only interested in the
>> superficial theory then that's enough.  The practical difficulty lies
>> in the fact that if you have noise in your system coupled with poor
>> numerical conditioning, it'll absolutely kill you.  In fact, even with
>> perfect measurements, doing the computation with fixed-precision math
>> will kill you, if the numerical conditioning is worse than the
>> available precision of your math.
> 
> This would somehow confirm my observation from above. So in theory,
> there is no fundamental lower bound.

If you model noise in your theory, then the fundamental lower bounds will 
make themselves evident.

Trying to draw conclusions from Matlab experimentation is a mistake.  
Using Matlab to point the way so you can sit down and grind through the 
math symbolically is a good idea, but you have to remember that Matlab is 
a _starting_ place, not an _ending_ place.

-- 
My liberal friends think I'm a conservative kook.
My conservative friends think I'm a liberal kook.
Why am I not happy that they have found common ground?

Tim Wescott, Communications, Control, Circuits & Software
http://www.wescottdesign.com

Hi,

Am 1/28/2013 8:04 AM, schrieb Vladimir Vassilevsky:
> "robert bristow-johnson" <rbj@audioimagination.com> wrote in message 
> news:ke66ml$dr3$1@dont-email.me...
>> On 1/28/13 10:27 AM, Vladimir Vassilevsky wrote:
>>> "Peter Mairhofer"<63832452@gmx.net>  wrote:
>>>
>>>> Is there a theoretical minimum for the sampling rate of the output y(t)
>>>> in a SysId setup? I know that the input signal x(t) must contain all
>>>> frequency components ("persistently exciting") of the unknown filter.
>>>
>>> Sample rate is not relevant.
>>>
>> well, the continuous-time system has "features" in the frequency response.
> 
> Any piece of frequency response contains full information about entire 
> frequency response.

Do you mean "any piece of system response"?

Can you detail this? "Any piece of of frequency reponse" is for me how
the amplitude and phase is changed by the system at one particular
frequency. This gives no information about the entire frequency range.

> [...]
> That is practical implementation issue.
> 
>>> Leaving aside misc. pathologies, you just need minimum of N samples to 
>>> solve
>>> for N coefficients.
>>
>> is this sysid model an FIR?

The actual system is an analog filter h(t); it is assumed that it can be
approximated well enough with a Q-tap FIR filter.

> [...]
>>> It doesn't matter how those samples are taken.
>>
>> i don't understand what variation in "how" you refer to, Vlad.
> 
> You have H(t) sampled at arbitrary time instants t1,t2....tN.
> That is enough to solve for N coefficients.

Ok, thanks.
That means taking the samples at the Nyquist rate somehow "ensures" that
the linear system is stable.

Dropping an arbitrary set of rows could yield a bad-conditioned matrix.

>> you might be able to ameliorate some aliasing problems in the frequency 
>> response.
> 
> No aliasing. If your model is adequate, the exact fit at t1,t2...tN means 
> exact fit everywhere.

As far as I understand (see my other response) aliasing is your friend
here (because it keeps the information). Low-Pass filtering before
(sub)-sampling would remove information.

Peter

Hi,

Thank you (and all) for the discussion.
I searched some typical books in adaptive filter theory but this
questions is never addressed.

Am 1/28/2013 9:20 AM, schrieb Tim Wescott:
> On Mon, 28 Jan 2013 00:09:51 -0800, Peter Mairhofer wrote:
> 
>> Hi,
>>
>> Is there a theoretical minimum for the sampling rate of the output y(t)
>> in a SysId setup?
> 
> The dangerous word in your question is "theoretical".  If you never need 
> a usable answer, then see Vladimir's response, because yes, in theory you 
> could do this at just about any sampling rate.

Ok, thank you. Acutally, I am interested in both cases - theoretical and
practical to understand my observations.

For the theoretical side, think about that: Is it a fundamental barrier
or could it beaten by a method which is not yet discovered? E.g., we
know now that the Nyquist rate is not the fundamental barrier for the
sampling rate, it's the information content (c.f. Compressive Sampling).
It is just a matter if it is possible to develop an appropriate sampling
scheme exploiting a specific structure.

So right now it sounds to me that there is no such fundamental limit
which has been proven. So it could be possible to discover completely
new ways of SysId, sampling the output at a lower rate which could even
be physically reliable.

My observation from simple MATLAB experiments is:

1) I can use any subset of the equations and get a solution. This is in
one sense obvious because it is LS. On the other hand, it contradicts my
intuition. Suppose my Nyquist rate for input x(t) and the unknown system
h(t) is 10kHz and h(t) operates (amongst others) also at 4 kHz. If I
sample the output y(t) at 5kHz, how should the reconstruction algorithm
know about the 4kHz range?

2) Using a lowpass filter to supress aliasing before sampling still
gives a valid result when the numbers are absolutely accurate. Minor
perturbations however, give absolutely wrong results.

My conclusion therefore is: The (implicit) aliasing when just using
*perfect* pointwise values of y(t) does the job. Low-pass filtering
removes the vital information from the system, rendering the
reconstruction very unstable.

(2) even gives correct results under perfect numerical conditions. My
mathematical interpretations for this is: The columns in the matrix
become nearly linearly dependent because of the "smoothing".

>> I know that the input signal x(t) must contain all
>> frequency components ("persistently exciting") of the unknown filter.
> 
> You know wrong. The input signal must contain _enough_ frequency 
> components to persistently excite the unknown filter. [...]

Sure. Sorry, that's what I meant.

>> When using a simple LS method to solve for the coefficients, the
>> overdetermined system of equations suggests that it is possible to drop
>> rows (as long as there are enough linearly independent) and thus
>> (implicitely) subsample y(t).
> 
> Yes, that is the suggestion.  And if you're only interested in the 
> superficial theory then that's enough.  The practical difficulty lies in 
> the fact that if you have noise in your system coupled with poor 
> numerical conditioning, it'll absolutely kill you.  In fact, even with 
> perfect measurements, doing the computation with fixed-precision math 
> will kill you, if the numerical conditioning is worse than the available 
> precision of your math.

This would somehow confirm my observation from above.
So in theory, there is no fundamental lower bound.

Thanks,
Peter

On Mon, 28 Jan 2013 00:09:51 -0800, Peter Mairhofer wrote:

> Hi,
> 
> Is there a theoretical minimum for the sampling rate of the output y(t)
> in a SysId setup?

The dangerous word in your question is "theoretical".  If you never need 
a usable answer, then see Vladimir's response, because yes, in theory you 
could do this at just about any sampling rate.

> I know that the input signal x(t) must contain all
> frequency components ("persistently exciting") of the unknown filter.

You know wrong.  The input signal must contain _enough_ frequency 
components to persistently excite the unknown filter.  For a system ID 
that's going to be carried out over a long time period, this essentially 
means that the input signal must contain as many frequency components as 
you have states in the filter.

> When using a simple LS method to solve for the coefficients, the
> overdetermined system of equations suggests that it is possible to drop
> rows (as long as there are enough linearly independent) and thus
> (implicitely) subsample y(t).

Yes, that is the suggestion.  And if you're only interested in the 
superficial theory then that's enough.  The practical difficulty lies in 
the fact that if you have noise in your system coupled with poor 
numerical conditioning, it'll absolutely kill you.  In fact, even with 
perfect measurements, doing the computation with fixed-precision math 
will kill you, if the numerical conditioning is worse than the available 
precision of your math.

(I just realized that I don't know how much work has been done on 
constructing LS system ID problems to avoid poor numerical conditioning!  
Stand by for a post...)

-- 
My liberal friends think I'm a conservative kook.
My conservative friends think I'm a liberal kook.
Why am I not happy that they have found common ground?

Tim Wescott, Communications, Control, Circuits & Software
http://www.wescottdesign.com

On Mon, 28 Jan 2013 10:04:31 -0600, Vladimir Vassilevsky wrote:

> "robert bristow-johnson" <rbj@audioimagination.com> wrote in message
> news:ke66ml$dr3$1@dont-email.me...
>> On 1/28/13 10:27 AM, Vladimir Vassilevsky wrote:
>>> "Peter Mairhofer"<63832452@gmx.net>  wrote:
>>>
>>>> Is there a theoretical minimum for the sampling rate of the output
>>>> y(t)
>>>> in a SysId setup? I know that the input signal x(t) must contain all
>>>> frequency components ("persistently exciting") of the unknown filter.
>>>
>>> Sample rate is not relevant.
>>>
>> well, the continuous-time system has "features" in the frequency
>> response.
> 
> Any piece of frequency response contains full information about entire
> frequency response.
> 
>>  i would hope that the sample rate is well over twice the frequency of
>>  the
>> highest-frequency feature.
> 
> That is practical implementation issue.

That's like saying that packing food for a long trip in the wilderness is 
a practical implementation issue.  Throw noise into the mix, and you're 
not going to _practically_ get good detail on high-frequency features of 
the system.

You could reasonably argue that if you're going to use the unknown system 
while sampling at F_s, that sampling at F_s will likely get you an 
adequate model, though.

-- 
Tim Wescott
Control system and signal processing consulting
www.wescottdesign.com

"robert bristow-johnson" <rbj@audioimagination.com> wrote in message 
news:ke66ml$dr3$1@dont-email.me...
> On 1/28/13 10:27 AM, Vladimir Vassilevsky wrote:
>> "Peter Mairhofer"<63832452@gmx.net>  wrote:
>>
>>> Is there a theoretical minimum for the sampling rate of the output y(t)
>>> in a SysId setup? I know that the input signal x(t) must contain all
>>> frequency components ("persistently exciting") of the unknown filter.
>>
>> Sample rate is not relevant.
>>
> well, the continuous-time system has "features" in the frequency response.

Any piece of frequency response contains full information about entire 
frequency response.

>  i would hope that the sample rate is well over twice the frequency of the 
> highest-frequency feature.

That is practical implementation issue.

>> Leaving aside misc. pathologies, you just need minimum of N samples to 
>> solve
>> for N coefficients.
>
> is this sysid model an FIR?

Does not matter as long as model is adequate.

>> It doesn't matter how those samples are taken.
>
> i don't understand what variation in "how" you refer to, Vlad.

You have H(t) sampled at arbitrary time instants t1,t2....tN.
That is enough to solve for N coefficients.

> you might be able to ameliorate some aliasing problems in the frequency 
> response.

No aliasing. If your model is adequate, the exact fit at t1,t2...tN means 
exact fit everywhere.

>> It won't
>> be practical but for classwork example, but that is a different question.
>
> i wasn't quite clear on the OP's application, but if it's a classwork 
> thing, and he has N contiguous input and output samples, i guess it's a 
> case of DFT, point-by-point division, and inverse DFT.

Nobody said the word "DFT" yet.


VLV

On 1/28/13 10:27 AM, Vladimir Vassilevsky wrote:
> "Peter Mairhofer"<63832452@gmx.net>  wrote:
>
>> Is there a theoretical minimum for the sampling rate of the output y(t)
>> in a SysId setup? I know that the input signal x(t) must contain all
>> frequency components ("persistently exciting") of the unknown filter.
>
> Sample rate is not relevant.
>

well, the continuous-time system has "features" in the frequency 
response.  i would hope that the sample rate is well over twice the 
frequency of the highest-frequency feature.

>> When using a simple LS method to solve for the coefficients, the
>> overdetermined system of equations suggests that it is possible to drop
>> rows (as long as there are enough linearly independent) and thus
>> (implicitely) subsample y(t).
>
> Leaving aside misc. pathologies, you just need minimum of N samples to solve
> for N coefficients.

is this sysid model an FIR?

> It doesn't matter how those samples are taken.

i don't understand what variation in "how" you refer to, Vlad.  if you 
decimate by regularly throwing away 9 out of every 10 samples, that is 
effectively cutting the sample rate by a factor of 10 and i would think 
that would matter quite a bit.  by having unequal spacing (but *known* 
spacing) of the samples, you might be able to ameliorate some aliasing 
problems in the frequency response.

> It won't
> be practical but for classwork example, but that is a different question.

i wasn't quite clear on the OP's application, but if it's a classwork 
thing, and he has N contiguous input and output samples, i guess it's a 
case of DFT, point-by-point division, and inverse DFT.

-- 

r b-j                  rbj@audioimagination.com

"Imagination is more important than knowledge."