DSPRelated.com
Forums

Decorrelation

Started by Tim Wescott August 31, 2016
I've got a nice real-world problem that I've reduced down to the 
following homework-y sounding problem:

Given a system h, I give it an excitation x in sampled time, and I get 
back

r(k) = h(x(k) + n1(k)) + n2(k)

where r(k) is a known received signal, x(k) is a presumed input signal, 
and n1(k) and n2(k) are unknown noise.  I'm taking a repeated signal for 
both x and r and averaging the snot out of them over the cycle, so the 
actual noise is presumed to be low.  However, the system isn't exactly 
linear, so n2(k) has some "noise" that's actually nonlinear effects.

I want to put on my "easy math" blinders and assume that h is linear and 
time invariant (it is, sorta), and that the noise signals are white and 
Gaussian.  Then I want to find the impulse response of h.

I'm currently doing this by absolute brute force.  I'm making a matrix of 
delayed versions of x and making the matrix

    [x1 0  0  ...]
    [x2 x1 0  ...]
X = [x3 x2 x1 ...]
    [     ...    ]

Then I'm calculating the least-mean-squared error solution to

r = X * b

Then I'm taking the vector b as the impulse response of h.

x and r are, at the moment, 4096 points long, and I'm getting 512 points 
for b.  So I'm calculating the above equation for a 4096 x 512 point 
matrix.

The way I'm doing it gives good results, but it seems exceedingly 
wasteful of computing resources.  It not only seems that this should be 
something that could be done handily using the FFT, but that I should 
know how to do it.  So far, I've failed.  It ought to be some variation 
of Wiener filtering -- but, as I've said, so far all I get when trying to 
apply FFTs to the problem is absolute crap coming back.

Is there a NAME for the FFT magic I'm trying to do, and is there a white 
paper (or a book chapter) on it that I can go read?  I'm not quite at the 
point of buying a chain saw and attempting to make a living by carving 
tree trunks into cheezy-looking bears and selling them to tourists, but 
I'm getting close.

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

I'm looking for work -- see my website!
On Wed, 31 Aug 2016 12:39:00 -0500, Tim Wescott
<seemywebsite@myfooter.really> wrote:

>I've got a nice real-world problem that I've reduced down to the >following homework-y sounding problem: > >Given a system h, I give it an excitation x in sampled time, and I get >back > >r(k) = h(x(k) + n1(k)) + n2(k) > >where r(k) is a known received signal, x(k) is a presumed input signal, >and n1(k) and n2(k) are unknown noise. I'm taking a repeated signal for >both x and r and averaging the snot out of them over the cycle, so the >actual noise is presumed to be low. However, the system isn't exactly >linear, so n2(k) has some "noise" that's actually nonlinear effects. > >I want to put on my "easy math" blinders and assume that h is linear and >time invariant (it is, sorta), and that the noise signals are white and >Gaussian. Then I want to find the impulse response of h. > >I'm currently doing this by absolute brute force. I'm making a matrix of >delayed versions of x and making the matrix > > [x1 0 0 ...] > [x2 x1 0 ...] >X = [x3 x2 x1 ...] > [ ... ] > >Then I'm calculating the least-mean-squared error solution to > >r = X * b > >Then I'm taking the vector b as the impulse response of h. > >x and r are, at the moment, 4096 points long, and I'm getting 512 points >for b. So I'm calculating the above equation for a 4096 x 512 point >matrix. > >The way I'm doing it gives good results, but it seems exceedingly >wasteful of computing resources. It not only seems that this should be >something that could be done handily using the FFT, but that I should >know how to do it. So far, I've failed. It ought to be some variation >of Wiener filtering -- but, as I've said, so far all I get when trying to >apply FFTs to the problem is absolute crap coming back. > >Is there a NAME for the FFT magic I'm trying to do, and is there a white >paper (or a book chapter) on it that I can go read? I'm not quite at the >point of buying a chain saw and attempting to make a living by carving >tree trunks into cheezy-looking bears and selling them to tourists, but >I'm getting close. > >-- > >Tim Wescott >Wescott Design Services >http://www.wescottdesign.com > >I'm looking for work -- see my website!
If you want a reasonable frequency response in place of the impulse response, or to compute the impulse response, then if the system cooperates you can sweep x(k) with a linear FM. If this is done coherently so that the phase of x(k) tracked across to r(k), then the FT of r(k) will essentially be the impulse response. The idea is to do the function of a Vector Signal Analyzer, where x(k) is the swept stimulus and r(k) is the measured output. This does assume that n1 and n2 are not screwing things up too badly, and that h(k) is reasonably LTI.
On Wed, 31 Aug 2016 18:07:45 GMT, eric.jacobsen@ieee.org (Eric
Jacobsen) wrote:

>On Wed, 31 Aug 2016 12:39:00 -0500, Tim Wescott ><seemywebsite@myfooter.really> wrote: > >>I've got a nice real-world problem that I've reduced down to the >>following homework-y sounding problem: >> >>Given a system h, I give it an excitation x in sampled time, and I get >>back >> >>r(k) = h(x(k) + n1(k)) + n2(k) >> >>where r(k) is a known received signal, x(k) is a presumed input signal, >>and n1(k) and n2(k) are unknown noise. I'm taking a repeated signal for >>both x and r and averaging the snot out of them over the cycle, so the >>actual noise is presumed to be low. However, the system isn't exactly >>linear, so n2(k) has some "noise" that's actually nonlinear effects. >> >>I want to put on my "easy math" blinders and assume that h is linear and >>time invariant (it is, sorta), and that the noise signals are white and >>Gaussian. Then I want to find the impulse response of h. >> >>I'm currently doing this by absolute brute force. I'm making a matrix of >>delayed versions of x and making the matrix >> >> [x1 0 0 ...] >> [x2 x1 0 ...] >>X = [x3 x2 x1 ...] >> [ ... ] >> >>Then I'm calculating the least-mean-squared error solution to >> >>r = X * b >> >>Then I'm taking the vector b as the impulse response of h. >> >>x and r are, at the moment, 4096 points long, and I'm getting 512 points >>for b. So I'm calculating the above equation for a 4096 x 512 point >>matrix. >> >>The way I'm doing it gives good results, but it seems exceedingly >>wasteful of computing resources. It not only seems that this should be >>something that could be done handily using the FFT, but that I should >>know how to do it. So far, I've failed. It ought to be some variation >>of Wiener filtering -- but, as I've said, so far all I get when trying to >>apply FFTs to the problem is absolute crap coming back. >> >>Is there a NAME for the FFT magic I'm trying to do, and is there a white >>paper (or a book chapter) on it that I can go read? I'm not quite at the >>point of buying a chain saw and attempting to make a living by carving >>tree trunks into cheezy-looking bears and selling them to tourists, but >>I'm getting close. >> >>-- >> >>Tim Wescott >>Wescott Design Services >>http://www.wescottdesign.com >> >>I'm looking for work -- see my website! > >If you want a reasonable frequency response in place of the impulse >response, or to compute the impulse response, then if the system >cooperates you can sweep x(k) with a linear FM. If this is done >coherently so that the phase of x(k) tracked across to r(k), then the >FT of r(k) will essentially be the impulse response. > >The idea is to do the function of a Vector Signal Analyzer, where x(k) >is the swept stimulus and r(k) is the measured output. This does >assume that n1 and n2 are not screwing things up too badly, and that >h(k) is reasonably LTI. > >
Looking at this again I may not have gotten your notation right. Arg. Regardless, IF you can get a coherent frequency response with a linear sweep (or somehow), then the impulse response is just the FT of that. I know you know that, though.
On Wed, 31 Aug 2016 18:07:45 +0000, Eric Jacobsen wrote:

> On Wed, 31 Aug 2016 12:39:00 -0500, Tim Wescott > <seemywebsite@myfooter.really> wrote: > >>I've got a nice real-world problem that I've reduced down to the >>following homework-y sounding problem: >> >>Given a system h, I give it an excitation x in sampled time, and I get >>back >> >>r(k) = h(x(k) + n1(k)) + n2(k) >> >>where r(k) is a known received signal, x(k) is a presumed input signal, >>and n1(k) and n2(k) are unknown noise. I'm taking a repeated signal for >>both x and r and averaging the snot out of them over the cycle, so the >>actual noise is presumed to be low. However, the system isn't exactly >>linear, so n2(k) has some "noise" that's actually nonlinear effects. >> >>I want to put on my "easy math" blinders and assume that h is linear and >>time invariant (it is, sorta), and that the noise signals are white and >>Gaussian. Then I want to find the impulse response of h. >> >>I'm currently doing this by absolute brute force. I'm making a matrix >>of delayed versions of x and making the matrix >> >> [x1 0 0 ...] >> [x2 x1 0 ...] >>X = [x3 x2 x1 ...] >> [ ... ] >> >>Then I'm calculating the least-mean-squared error solution to >> >>r = X * b >> >>Then I'm taking the vector b as the impulse response of h. >> >>x and r are, at the moment, 4096 points long, and I'm getting 512 points >>for b. So I'm calculating the above equation for a 4096 x 512 point >>matrix. >> >>The way I'm doing it gives good results, but it seems exceedingly >>wasteful of computing resources. It not only seems that this should be >>something that could be done handily using the FFT, but that I should >>know how to do it. So far, I've failed. It ought to be some variation >>of Wiener filtering -- but, as I've said, so far all I get when trying >>to apply FFTs to the problem is absolute crap coming back. >> >>Is there a NAME for the FFT magic I'm trying to do, and is there a white >>paper (or a book chapter) on it that I can go read? I'm not quite at >>the point of buying a chain saw and attempting to make a living by >>carving tree trunks into cheezy-looking bears and selling them to >>tourists, but I'm getting close. >> >>-- >> >>Tim Wescott Wescott Design Services http://www.wescottdesign.com >> >>I'm looking for work -- see my website! > > If you want a reasonable frequency response in place of the impulse > response, or to compute the impulse response, then if the system > cooperates you can sweep x(k) with a linear FM. If this is done > coherently so that the phase of x(k) tracked across to r(k), then the FT > of r(k) will essentially be the impulse response. > > The idea is to do the function of a Vector Signal Analyzer, where x(k) > is the swept stimulus and r(k) is the measured output. This does > assume that n1 and n2 are not screwing things up too badly, and that > h(k) is reasonably LTI.
I don't get to choose x(k). I've got a huge collection of raw signal recordings from a customer, from which I'm trying to reconstruct channel models. I suppose I should just do it by brute force and be happy... -- Tim Wescott Control systems, embedded software and circuit design I'm looking for work! See my website if you're interested http://www.wescottdesign.com
On Wed, 31 Aug 2016 15:51:04 -0500, Tim Wescott <tim@seemywebsite.com>
wrote:

>On Wed, 31 Aug 2016 18:07:45 +0000, Eric Jacobsen wrote: > >> On Wed, 31 Aug 2016 12:39:00 -0500, Tim Wescott >> <seemywebsite@myfooter.really> wrote: >> >>>I've got a nice real-world problem that I've reduced down to the >>>following homework-y sounding problem: >>> >>>Given a system h, I give it an excitation x in sampled time, and I get >>>back >>> >>>r(k) = h(x(k) + n1(k)) + n2(k) >>> >>>where r(k) is a known received signal, x(k) is a presumed input signal, >>>and n1(k) and n2(k) are unknown noise. I'm taking a repeated signal for >>>both x and r and averaging the snot out of them over the cycle, so the >>>actual noise is presumed to be low. However, the system isn't exactly >>>linear, so n2(k) has some "noise" that's actually nonlinear effects. >>> >>>I want to put on my "easy math" blinders and assume that h is linear and >>>time invariant (it is, sorta), and that the noise signals are white and >>>Gaussian. Then I want to find the impulse response of h. >>> >>>I'm currently doing this by absolute brute force. I'm making a matrix >>>of delayed versions of x and making the matrix >>> >>> [x1 0 0 ...] >>> [x2 x1 0 ...] >>>X = [x3 x2 x1 ...] >>> [ ... ] >>> >>>Then I'm calculating the least-mean-squared error solution to >>> >>>r = X * b >>> >>>Then I'm taking the vector b as the impulse response of h. >>> >>>x and r are, at the moment, 4096 points long, and I'm getting 512 points >>>for b. So I'm calculating the above equation for a 4096 x 512 point >>>matrix. >>> >>>The way I'm doing it gives good results, but it seems exceedingly >>>wasteful of computing resources. It not only seems that this should be >>>something that could be done handily using the FFT, but that I should >>>know how to do it. So far, I've failed. It ought to be some variation >>>of Wiener filtering -- but, as I've said, so far all I get when trying >>>to apply FFTs to the problem is absolute crap coming back. >>> >>>Is there a NAME for the FFT magic I'm trying to do, and is there a white >>>paper (or a book chapter) on it that I can go read? I'm not quite at >>>the point of buying a chain saw and attempting to make a living by >>>carving tree trunks into cheezy-looking bears and selling them to >>>tourists, but I'm getting close. >>> >>>-- >>> >>>Tim Wescott Wescott Design Services http://www.wescottdesign.com >>> >>>I'm looking for work -- see my website! >> >> If you want a reasonable frequency response in place of the impulse >> response, or to compute the impulse response, then if the system >> cooperates you can sweep x(k) with a linear FM. If this is done >> coherently so that the phase of x(k) tracked across to r(k), then the FT >> of r(k) will essentially be the impulse response. >> >> The idea is to do the function of a Vector Signal Analyzer, where x(k) >> is the swept stimulus and r(k) is the measured output. This does >> assume that n1 and n2 are not screwing things up too badly, and that >> h(k) is reasonably LTI. > >I don't get to choose x(k). I've got a huge collection of raw signal >recordings from a customer, from which I'm trying to reconstruct channel >models. I suppose I should just do it by brute force and be happy...
Ouch. Can't do any better than whatever the characteristics of the existing x(k)s were then, so the results could be sparse. In the past I've done something similar to this: as long as the LTI assumption holds reasonably, you can everage the crap out of the spectral content of r(k) to get an idea of the flatness or lack thereof of the frequency response. A rough estimate of an impulse response can then be made, but without a means to keep coherence between all of the samples it gets pretty hand-wavy. Sometimes that's as good as you can do.
Tim Wescott wrote:
> I've got a nice real-world problem that I've reduced down to the > following homework-y sounding problem: > > Given a system h, I give it an excitation x in sampled time, and I get > back > > r(k) = h(x(k) + n1(k)) + n2(k) > > where r(k) is a known received signal, x(k) is a presumed input signal, > and n1(k) and n2(k) are unknown noise. I'm taking a repeated signal for > both x and r and averaging the snot out of them over the cycle, so the > actual noise is presumed to be low. However, the system isn't exactly > linear, so n2(k) has some "noise" that's actually nonlinear effects. > > I want to put on my "easy math" blinders and assume that h is linear and > time invariant (it is, sorta), and that the noise signals are white and > Gaussian. Then I want to find the impulse response of h. > > I'm currently doing this by absolute brute force. I'm making a matrix of > delayed versions of x and making the matrix > > [x1 0 0 ...] > [x2 x1 0 ...] > X = [x3 x2 x1 ...] > [ ... ] > > Then I'm calculating the least-mean-squared error solution to > > r = X * b > > Then I'm taking the vector b as the impulse response of h. > > x and r are, at the moment, 4096 points long, and I'm getting 512 points > for b. So I'm calculating the above equation for a 4096 x 512 point > matrix. > > The way I'm doing it gives good results, but it seems exceedingly > wasteful of computing resources. It not only seems that this should be > something that could be done handily using the FFT, but that I should > know how to do it. So far, I've failed. It ought to be some variation > of Wiener filtering -- but, as I've said, so far all I get when trying to > apply FFTs to the problem is absolute crap coming back. > > Is there a NAME for the FFT magic I'm trying to do, and is there a white > paper (or a book chapter) on it that I can go read? I'm not quite at the > point of buying a chain saw and attempting to make a living by carving > tree trunks into cheezy-looking bears and selling them to tourists, but > I'm getting close. >
Is there a method involving deconvolution lurking under all that? -- Les Cargill
Am Mittwoch, 31. August 2016 19:39:07 UTC+2 schrieb Tim Wescott:
> > The way I'm doing it gives good results, but it seems exceedingly > wasteful of computing resources. It not only seems that this should be > something that could be done handily using the FFT, but that I should > know how to do it. So far, I've failed. It ought to be some variation > of Wiener filtering -- but, as I've said, so far all I get when trying to > apply FFTs to the problem is absolute crap coming back.
Did you try Pwelch's method for estimating the transfer function? (See Matlab's tfestimate. It's FFT-based). It's fast and robust. But I suppose you might need more than 4096 samples if you expect the impulse response to span about 512 samples. You might need much more to let the method work on large blocks and have enough blocks for the averaging. Another approach would be to window and zero-pad x and b so that the columns of X contain the same non-zero coefficients only shifted. Then you could try to solve the normal equation system instead: A r = X^T b with A = X^T X Due to the structure of X, A will be a symmetric Toeplitz matrix (and hopefully a positive definite one). For these kinds of matrices there are special O(n^2) solvers. Computing A can be done using an FFT since it contains the scaled auto-correlation of x. And X^T b is a cross- correlation between x and b which can also be computed using the FFT. Though, numerically, this approach is worse as it basically involves squaring X. So, if the condition number of X is already bad, this is probably a bad idea. But it'll be faster. The condition number of A is basically close to the ratio of the highest to the lowest power spectral density of x. If your signal x is white or whitish then the condition number of A should be fine and you shouldn't have any numerical problems. Solving similar Toeplitz systems is what many speech coding people do for estimating the all-pole filter that turns a buzzing sound or noise into speech. Cheers! SG
On 09/01/2016 04:51 AM, Tim Wescott wrote:
> On Wed, 31 Aug 2016 18:07:45 +0000, Eric Jacobsen wrote: > >> On Wed, 31 Aug 2016 12:39:00 -0500, Tim Wescott >> <seemywebsite@myfooter.really> wrote: >> >>> I've got a nice real-world problem that I've reduced down to the >>> following homework-y sounding problem: >>> >>> Given a system h, I give it an excitation x in sampled time, and I get >>> back >>> >>> r(k) = h(x(k) + n1(k)) + n2(k) >>> >>> where r(k) is a known received signal, x(k) is a presumed input signal, >>> and n1(k) and n2(k) are unknown noise. I'm taking a repeated signal for >>> both x and r and averaging the snot out of them over the cycle, so the >>> actual noise is presumed to be low. However, the system isn't exactly >>> linear, so n2(k) has some "noise" that's actually nonlinear effects. >>> >>> I want to put on my "easy math" blinders and assume that h is linear and >>> time invariant (it is, sorta), and that the noise signals are white and >>> Gaussian. Then I want to find the impulse response of h. >>> >>> I'm currently doing this by absolute brute force. I'm making a matrix >>> of delayed versions of x and making the matrix >>> >>> [x1 0 0 ...] >>> [x2 x1 0 ...] >>> X = [x3 x2 x1 ...] >>> [ ... ] >>> >>> Then I'm calculating the least-mean-squared error solution to >>> >>> r = X * b >>> >>> Then I'm taking the vector b as the impulse response of h. >>> >>> x and r are, at the moment, 4096 points long, and I'm getting 512 points >>> for b. So I'm calculating the above equation for a 4096 x 512 point >>> matrix. >>> >>> The way I'm doing it gives good results, but it seems exceedingly >>> wasteful of computing resources. It not only seems that this should be >>> something that could be done handily using the FFT, but that I should >>> know how to do it. So far, I've failed. It ought to be some variation >>> of Wiener filtering -- but, as I've said, so far all I get when trying >>> to apply FFTs to the problem is absolute crap coming back. >>> >>> Is there a NAME for the FFT magic I'm trying to do, and is there a white >>> paper (or a book chapter) on it that I can go read? I'm not quite at >>> the point of buying a chain saw and attempting to make a living by >>> carving tree trunks into cheezy-looking bears and selling them to >>> tourists, but I'm getting close. >>> >>> -- >>> >>> Tim Wescott Wescott Design Services http://www.wescottdesign.com >>> >>> I'm looking for work -- see my website! >> >> If you want a reasonable frequency response in place of the impulse >> response, or to compute the impulse response, then if the system >> cooperates you can sweep x(k) with a linear FM. If this is done >> coherently so that the phase of x(k) tracked across to r(k), then the FT >> of r(k) will essentially be the impulse response. >> >> The idea is to do the function of a Vector Signal Analyzer, where x(k) >> is the swept stimulus and r(k) is the measured output. This does >> assume that n1 and n2 are not screwing things up too badly, and that >> h(k) is reasonably LTI. > > I don't get to choose x(k). I've got a huge collection of raw signal > recordings from a customer, from which I'm trying to reconstruct channel > models. I suppose I should just do it by brute force and be happy... >
How broadband are those stimuli, either individually or when you look at their sum? That's usually what determines how well you can characterise the system. Regards, Steve
On Thursday, September 1, 2016 at 7:52:22 AM UTC+2, SG wrote:
> [...] > Another approach would be to window and zero-pad x and b so that the > columns of X contain the same non-zero coefficients only shifted. Then > you could try to solve the normal equation system instead: > > A r = X^T b with A = X^T X > > Due to the structure of X, A will be a symmetric Toeplitz matrix (and > hopefully a positive definite one). For these kinds of matrices there > are special O(n^2) solvers. Computing A can be done using an FFT since > it contains the scaled auto-correlation of x. And X^T b is a cross- > correlation between x and b which can also be computed using the FFT. > > Though, numerically, this approach is worse as it basically involves > squaring X. So, if the condition number of X is already bad, this is > probably a bad idea. But it'll be faster. The condition number of A is > basically close to the ratio of the highest to the lowest power > spectral density of x. If your signal x is white or whitish then the > condition number of A should be fine and you shouldn't have any > numerical problems. Solving similar Toeplitz systems is what many > speech coding people do for estimating the all-pole filter that turns > a buzzing sound or noise into speech.
Warning: I might have gotton your b and r mixed up. I was assuming r to be your impulse response and b to be the convolved signal: conv(x + n1, r) + n2 = b And there might be a problem with not knowing n1. Your approach and the ones I was talking about basically ignore n1 completely. This may be a bad idea. In the spectral domain with k overlapping windows you should be able to treat all the frequencies in isolation. So, for some frequency you collect all the respective FFT bins into columns vectors x, n1, n2, b. Then you have the following overdetermined equation system for each frequency bin: (x + n1) r + n2 = b (1) where r is a complex scalar (unknown amplitude and phase response). These are k equations (due to k FFT blocks) with a single unknown r. Solving this in a least squares sense can be done by the corresponding normal equation: (x + n1)^H (x + n1) r = (x + n1)^H (b - n2) (2) Assuming n2 does not correlate with x nor with n1 we can ignore it: (x + n1)^H (x + n1) r = (x + n1)^H b (3) Assuming x does not correlate with n1 we can rewrite this into (x^H x + n1^H n1) r = (x + n1)^H b (4) This is a little nicer. n1^H n1 is basically proportional to the power spectral density of n1 which you might know or be able to estimate somehow. But this still leaves n1 on the right hand side. I don't see a an easy way to get rid of it -- unless you're fine with solving a nonlinear equation system that tries to recover n1 as well constrained with n1^H n1 = some_power_level.
> Cheers! > SG
Am Donnerstag, 1. September 2016 10:58:13 UTC+2 schrieb SG:
> [...] > In the spectral domain with k overlapping windows you should be able > to treat all the frequencies in isolation. So, for some frequency you > collect all the respective FFT bins into columns vectors x, n1, n2, b. > Then you have the following overdetermined equation system for each > frequency bin: > > (x + n1) r + n2 = b (1) > > where r is a complex scalar (unknown amplitude and phase response). > These are k equations (due to k FFT blocks) with a single unknown r. > Solving this in a least squares sense can be done by the corresponding > normal equation: > > (x + n1)^H (x + n1) r = (x + n1)^H (b - n2) (2) > > Assuming n2 does not correlate with x nor with n1 we can ignore it: > > (x + n1)^H (x + n1) r = (x + n1)^H b (3) > > Assuming x does not correlate with n1 we can rewrite this into > > (x^H x + n1^H n1) r = (x + n1)^H b (4) > > This is a little nicer. n1^H n1 is basically proportional to the power > spectral density of n1 which you might know or be able to estimate > somehow. But this still leaves n1 on the right hand side. I don't see > a an easy way to get rid of it -- unless you're fine with solving a > nonlinear equation system that tries to recover n1 as well constrained > with n1^H n1 = some_power_level.
On second thought, this isn't nonlinear but underdetermined. You could select any n1 constrained to n1^H n1 = some_power_level which gives you (x + n1)^H b r = ------------------------- (x^H x + some_power_level) This tells you how well r is determined since you don't know n1 but only its norm, sqrt(some_power_level). And restoring your use of variable names: (x + n1)^H r h = ------------------------- (x^H x + some_power_level) Sorry for messing up r, h, b before. Your equation system confused me to think b is the received signal and r is the impulse response.
> > Cheers! > > SG