On Thu, 21 Dec 2006 17:44:00 -0500, Randy Yates <yates@ieee.org> wrote:
>What I'm wondering is if it can be proven that the algorithm in
>general will produce a stable IIR transfer function, i.e., how can you
>guarantee the poles it generates are always inside the unit circle?
As far as I know, you can't. However, in my experience, with only some
pathological exceptions, stable prototype systems have always yielded
stable models and unstable prototype systems have always yielded
unstable models. The exceptions were cases where pole(s) were extremely
close to the jw-axis (like real part > -0.0001) or to the unit circle
(like radius > .9999), where the model pole might slip to an unstable
position. That's why I was so surprised when the A-weighting filter
produced an unstable model -- it's a fairly benign transfer function.
The problem in this case had to do with the way that I implemented the
"artificial delay" mentioned in the magazine article. While the
explanation is probably too complex for a brief comp.dsp post, I can
summarize:
I use a "cosine" input instead of a "sine" input because of
observability problems with "sine" near half the sampling frequency.
(This is also mentioned in the article.) If the frequency of a sine
wave is "Fs/2", and you sample it at "Fs", you only capture the
zero-crossings. So, assuming that the output of the system at Fs/2 is
nonzero, the algorithm interprets the situation as
"zero-input/nonzero-output". It deals with this by placing one or more
poles at an angle of PI on the Z-plane, and those poles are often
outside the unit circle.
In my computer program I made the mistake of adding the "artificial
delay" not only to the output sequence, but to the input sequence as
well. Well, if you delay a cosine wave you phase-shift it. Trig
identities will tell you that a phase-shifted cosine is equal to a
cosine plus a sine. That little bit of sine wave on the input led to
the unstable pole(s) in the transfer function. I solved the problem by
eliminating the phase shift from the input sequence, allowing it only in
the output sequence.
The net effect of this is that if you inject "delta" samples of
"artificial delay" into the frequency response, the first "delta"
coefficients of the numerator of the transfer function will be
approximately zero. (Approximately, not exactly, because of the least
squares nature of the model fit.)
Now a teaser, for future study: It can be shown that a cosine input
formulation models the real part of the frequency response, and that a
sine input formulation models the imaginary part (and is almost always
unstable because of the observability problems mentioned above). When I
was in graduate school I was trying to find a way to eliminate that sine
formulation instability. But I had to move on with my life before I
solved the problem.
Greg
Reply by Randy Yates●December 21, 20062006-12-21
Greg Berchin <gberchin@comicast.net> writes:
> On 21 Dec 2006 14:10:13 -0800, "robert bristow-johnson"
> <rbj@audioimagination.com> wrote:
>
>>so what did you have to do? reflect the unstable poles into the
>>interior of the unit circle?
>
> Naw, I just implemented the algorithm CORRECTLY.
Hi Greg,
What I hear you saying is that, for a specific design, you had
a problem in your implementation that resulted in an unstable
set of coefficients.
What I'm wondering is if it can be proven that the algorithm in
general will produce a stable IIR transfer function, i.e., how can you
guarantee the poles it generates are always inside the unit circle?
--
% Randy Yates % "And all that I can do
%% Fuquay-Varina, NC % is say I'm sorry,
%%% 919-577-9882 % that's the way it goes..."
%%%% <yates@ieee.org> % Getting To The Point', *Balance of Power*, ELO
http://home.earthlink.net/~yatescr
Reply by Greg Berchin●December 21, 20062006-12-21
On 21 Dec 2006 14:10:13 -0800, "robert bristow-johnson"
<rbj@audioimagination.com> wrote:
>so what did you have to do? reflect the unstable poles into the
>interior of the unit circle?
Naw, I just implemented the algorithm CORRECTLY.
Greg
Reply by Greg Berchin●December 21, 20062006-12-21
As promised, here's a STABLE 48kHz set. See other thread if interested
in explanation as to why the previous set was unstable.
This set is designed for best performance between 20 Hz and 20 kHz.
-- Greg
B =
0.01147155239724
-0.0248548824653166
0.0323052801494283
-0.0555571372813964
0.233784933167266
0.392882400411297
-0.633249911806281
-0.479494344932334
0.322061493940389
0.197659563091056
0.00299879461389451
A =
1
-0.925699454182466
-0.992471193943543
0.837650096562845
0.22307303912603
-0.158404327757755
0.0184103295763937
Reply by robert bristow-johnson●December 21, 20062006-12-21
Greg Berchin wrote:
> As promised, here's a STABLE 48kHz set. See other thread if interested
> in explanation as to why the previous set was unstable.
>
> This set is designed for best performance between 20 Hz and 20 kHz.
so what did you have to do? reflect the unstable poles into the
interior of the unit circle?
r b-j
Reply by ●December 19, 20062006-12-19
robert bristow-johnson <rbj@audioimagination.com> wrote:
>
> Gerold Schrutz wrote:
>> Hello Bob,
>>
>> I'll send it to you via email.
>>
>
> why not post it? c'mon! spread the joy. it's the christmas season.
On 18 Dec 2006 04:05:04 -0800, "Robert Adams" <robert.adams@analog.com>
wrote:
>Presumably the order of this filter is low enough that one of the
>factoring programs would work on it (to get cascaded biquads);
That's actually how I discovered that it was unstable -- I did a partial
fraction decomposition and found to my horror that some of the poles lie
outside the unit circle.
Greg
Reply by Robert Adams●December 18, 20062006-12-18
Presumably the order of this filter is low enough that one of the
factoring programs would work on it (to get cascaded biquads); if so,
then I would guess that some of the large coefficients causing rbj's
incontinence issues would be dramatically lower.
Bob
Greg Berchin wrote:
Reply by Gerold Schrutz●December 18, 20062006-12-18
"robert bristow-johnson" <rbj@audioimagination.com> schrieb im Newsbeitrag
news:1166377947.282558.324210@t46g2000cwa.googlegroups.com...
>
> Gerold Schrutz wrote:
>> Hello Bob,
>>
>> I'll send it to you via email.
>>
>
> why not post it? c'mon! spread the joy. it's the christmas season.
>
> r b-j
>
Hello Robert,
I've got an C++ Object with coefs, filter structure, ..., so its hard to
post it!
If you want to have it, just give me a note.
Regards
Gerold
Reply by Philip Martel●December 17, 20062006-12-17
"Jerry Avins" <jya@ieee.org> wrote in message
news:9r-dnZb60rzTthnYnZ2dnUVZ_u2mnZ2d@rcn.net...
> Al Clark wrote:
>> Jerry Avins <jya@ieee.org> wrote in
>> news:RoWdnS_wtccSihnYnZ2dnUVZ_qDinZ2d@rcn.net:
>>
>>> Greg Berchin wrote:
>>>> On 15 Dec 2006 22:23:53 -0800, "robert bristow-johnson"
>>>> <rbj@audioimagination.com> wrote:
>>>>
>>>>> it oughta be okay to just apply the BLT* to the poles/zeros to get a
>>>>> digital filter.
>>>> Depends upon the sampling rate. That 12200 Hz pole pair makes life
>>>> very
>>>> difficult for the Bilinear Transform method if the sampling is 44100.
>>> What about upsampling before the filter and downsampling after it?
>>> Anyhow, is there a phase spec too? I see no reason for one in a
>>> weighting scheme.
>>>
>>> Jerry
>>
>> This has always been one of my solution approaches. It costs a bit more
>> in computation. You still need to over sample by quite a bit if you want
>> the curve to be very accurate since warping is still significant at
>> 20-30KHz.
>
> Prewarping will put the corner frequencies in the right place. The problem
> lies in the curve between them. It probably won't depart too severely even
> with a 3:1 oversampling rate, and 2:1 might be good enough. Has anyone
> tried? Experience beats speculation every time.
>
I would speculate that you are right, but...
> Jerry
> --
> Engineering is the art of making what you want from things you can get.
> �����������������������������������������������������������������������