>Now as to the "artificial delay"; when formulated using cosine waves,
>as in the magazine article and the Matlab code, the algorithm only
>models the real part of the frequency response.  Similarly, if
>formulated using sine waves instead of cosine waves, it only models
>the imaginary part.  (Unfortunately, a sine formulation has some
>stability problems related to the fact that sine waves whose frequency
>is close to half the sampling frequency are sampled very close to the
>zero crossings, so creating a sine formulation is a nontrivial task.)
>The easiest way to make some of the imaginary part of the frequency
>response real (and some of the real part imaginary) is to add a small
>amount of delay to the transfer function.  That's what the "artificial
>delay" parameter is for.  Many times you'll get a good fit without it.
>Other times one or two samples of delay will give you a better fit. If
>you need more than a few samples of delay, then something else is
>probably wrong.
>
Ah, this explains a great deal, Greg. Thank you!

- Terry

On Tue, 21 Jul 2009 09:45:23 -0500, "montlick"
<terry@tmlaboratories.com> wrote:

>Isn't this the same as setting all the A_m coefficients to 1, which can be
>done without modifying the program? Or am I not understanding something
>fundamental here? I have already tried this, and got no phase match and
>unrelated numerator and denominator coefficients.

I didn't quite follow Andor's explanation because of nomenclature
problems in this fixed text format.  But basically, you should be able
to start with the fundamental FDLS transfer function,

                      -1            -n
               b + b z   + ... + b z
        y(z)    0   1             n
        ---- = -----------------------
        u(z)          -1            -d
                1 + a z  + ... + a z
                     1            d

and CONSTRAIN the coefficients to be of the form you desire (numerator
coefficients the same as denominator coefficients, but in reverse
order):

                        -1          -n
               a + a   z   + ... + z
        y(z)    d   d-1
        ---- = -----------------------
        u(z)          -1            -d
                1 + a z  + ... + a z
                     1            d

Then just follow the rest of the procedure exactly as described in the
article, except instead of there being "a" and "b" coefficients, there
will just be "a" coefficients used in both the denominator and the
numerator.

I suspect that this is what Andor was trying to say.

>Also, I am confused about the proper usage of the argument "artificial
>delay" (the variable named "delta") in this program. I used the M=8 FDLS
>Example written up by Greg Berchin and Richard Lyons. I found that the
>accuracy of the result was very dependent on the value of delta. A value of
>0.001 seemed to give reasonably good accuracy. But other values did not.
>How does one choose the correct value a priori?

First, make certain that you have version 2.00 of the Matlab code.
Earlier versions implemented the delay in a slightly different manner
that can lead to unstable models being created from stable transfer
functions.

Now as to the "artificial delay"; when formulated using cosine waves,
as in the magazine article and the Matlab code, the algorithm only
models the real part of the frequency response.  Similarly, if
formulated using sine waves instead of cosine waves, it only models
the imaginary part.  (Unfortunately, a sine formulation has some
stability problems related to the fact that sine waves whose frequency
is close to half the sampling frequency are sampled very close to the
zero crossings, so creating a sine formulation is a nontrivial task.)
The easiest way to make some of the imaginary part of the frequency
response real (and some of the real part imaginary) is to add a small
amount of delay to the transfer function.  That's what the "artificial
delay" parameter is for.  Many times you'll get a good fit without it.
Other times one or two samples of delay will give you a better fit. If
you need more than a few samples of delay, then something else is
probably wrong.

By the way, the value of "delta" was intended to be an integer (a
value of "1" representing a delay of one sample period, for example),
but now that you mention it there is no mathematical reason that it
HAS to be.

Greg

Thank you for your reply, Andor.

>On 20 Jul., 14:06, Terry wrote:
>> Hi Folks,
>>
>> I came across this thread while searching for an alternative to the
Matlab
>> function iirgrpdelay. I don't have MATLAB, but use Octave instead.
>> Unfortunately, iirgrpdelay uses a compiled MATLAB MEX file for its
guts.
>>
>> I tried Greg's cool function (which is 100% Octave compatible!), but
it
>> doesn't give me quite what I need. I need a true allpass filter
(numerator
>> coefficients same as reversed denominator coefficients), like
iirgrpdelay
>> produces. Putting in all 1's for magnitude does not yield this.
>
>Hi Terry
>
>It is quite easy to modify Greg's FDLS to design proper allpass
>filters. Remember that the underlying equations are of the form
>
>A_m cos(w_m + phi_m) = - A_m sum_{k=1}^p a_k cos(-w_m*k + phi_m) + sum_
>{k=0}^z b_k cos(-w_m*k)
>
>where w_m, m=1,2,...M, are the frequencies where you sample the
>allpass frequncy response, A_m the amplitude and phi_m the phase at
>that frequency (p is the order of the recursive part, z the order of
>the transversal part). For an allpass filter of order p, the above
>equation becomes
>
>cos(w_m + phi_m) = -sum_{k=1}^p a_k cos(-w_m*k + phi_m) + sum_{k=0}^p
>a_{p-k} cos(-w_m*k),
>
>so you get
>
>cos(w_m + phi_m) - cos(-w_m*p) = sum_{k=1}^p a_k ( cos(-w_m*k + phi_m)
>+ cos(-w_k*(p-k)) )
>
>and now solve for the allpass coefficients a_k.
>
>Regards,
>Andor

Isn't this the same as setting all the A_m coefficients to 1, which can be
done without modifying the program? Or am I not understanding something
fundamental here? I have already tried this, and got no phase match and
unrelated numerator and denominator coefficients.

Also, I am confused about the proper usage of the argument "artificial
delay" (the variable named "delta") in this program. I used the M=8 FDLS
Example written up by Greg Berchin and Richard Lyons. I found that the
accuracy of the result was very dependent on the value of delta. A value of
0.001 seemed to give reasonably good accuracy. But other values did not.
How does one choose the correct value a priori?

Thanks,
Terry

On 20 Jul., 14:06, Terry wrote:
> Hi Folks,
>
> I came across this thread while searching for an alternative to the Matlab
> function iirgrpdelay. I don't have MATLAB, but use Octave instead.
> Unfortunately, iirgrpdelay uses a compiled MATLAB MEX file for its guts.
>
> I tried Greg's cool function (which is 100% Octave compatible!), but it
> doesn't give me quite what I need. I need a true allpass filter (numerator
> coefficients same as reversed denominator coefficients), like iirgrpdelay
> produces. Putting in all 1's for magnitude does not yield this.

Hi Terry

It is quite easy to modify Greg's FDLS to design proper allpass
filters. Remember that the underlying equations are of the form

A_m cos(w_m + phi_m) = - A_m sum_{k=1}^p a_k cos(-w_m*k + phi_m) + sum_
{k=0}^z b_k cos(-w_m*k)

where w_m, m=1,2,...M, are the frequencies where you sample the
allpass frequncy response, A_m the amplitude and phi_m the phase at
that frequency (p is the order of the recursive part, z the order of
the transversal part). For an allpass filter of order p, the above
equation becomes

cos(w_m + phi_m) = -sum_{k=1}^p a_k cos(-w_m*k + phi_m) + sum_{k=0}^p
a_{p-k} cos(-w_m*k),

so you get

cos(w_m + phi_m) - cos(-w_m*p) = sum_{k=1}^p a_k ( cos(-w_m*k + phi_m)
+ cos(-w_k*(p-k)) )

and now solve for the allpass coefficients a_k.

Regards,
Andor

>
> I've already attempted to use Markus Lang's old MATLAB program,
> ap_lgs_z.m, whose method is decribed in his IEEE papers "Allpass Filter
> Design and Applications" and "Simple and Robust Method for the Design of
> Allpass Filters Using Least-Squares Phase Error Criterion." But this
> required some changes to get it to run at all (his variable names are in
> German, so I had to do some guessing!). While it converges, I don't seem to
> get the right filter coefficients. And Marcus is working back in Germany
> and on one of those extended summer vacations that Europeans get, so I
> can't contact him.

Don't know about those extended summer vacations :-). But then,
Switzerland aint Europe.

Hi Folks,

I came across this thread while searching for an alternative to the Matlab
function iirgrpdelay. I don't have MATLAB, but use Octave instead.
Unfortunately, iirgrpdelay uses a compiled MATLAB MEX file for its guts.

I tried Greg's cool function (which is 100% Octave compatible!), but it
doesn't give me quite what I need. I need a true allpass filter (numerator
coefficients same as reversed denominator coefficients), like iirgrpdelay
produces. Putting in all 1's for magnitude does not yield this.

I've already attempted to use Markus Lang's old MATLAB program,
ap_lgs_z.m, whose method is decribed in his IEEE papers "Allpass Filter
Design and Applications" and "Simple and Robust Method for the Design of
Allpass Filters Using Least-Squares Phase Error Criterion." But this
required some changes to get it to run at all (his variable names are in
German, so I had to do some guessing!). While it converges, I don't seem to
get the right filter coefficients. And Marcus is working back in Germany
and on one of those extended summer vacations that Europeans get, so I
can't contact him.

Any help on a finding working all-pass filter design code would be greatly
appreciated! If it is in pretty much any other computer language and
well-commented (with examples), I can translate it to MATLAB/Octave.

Thanks,
Terry Montlick
(acoustics and occasional signal processing stuff)

On Sat, 19 Jan 2008 13:28:56 -0500, Jerry Avins <jya@ieee.org> wrote:

>Watch out! Some yahoo in the gummint will see it as code for 
>"commiecast" and come looking for you.

"When everyone is out to get you, paranoid is just good thinking." 
	~Dr. Johnny Fever

Greg Berchin wrote:
> On Tue, 15 Jan 2008 18:57:38 GMT, Al Clark <aclark@danvillesignal.com> wrote:
> 
>> If anyone else has resources to share, we would be glad to post it since 
>> the old dspguru has been static for so long.
> 
> I have a "C+" version of FDLS working.  I call it "C+" because it was
> originally written in Fortran, then converted to C, and now it's somewhere
> between C and C++.  You'll need a C++ compiler.  I tried to structure it as
> similarly as I could to the Matlab 2.00 version.  Because it uses a different
> method to compute the pseudoinverse, it produces results that are comparable
> to, but not identical to, the Matlab version.
> 
> Any volunteers for beta testing?  Contact me by removing the "i" from
> "comicast".

Watch out! Some yahoo in the gummint will see it as code for 
"commiecast" and come looking for you.

Jerry
-- 
Engineering is the art of making what you want from things you can get.
&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;&#4294967295;

On Tue, 15 Jan 2008 18:57:38 GMT, Al Clark <aclark@danvillesignal.com> wrote:

>If anyone else has resources to share, we would be glad to post it since 
>the old dspguru has been static for so long.

I have a "C+" version of FDLS working.  I call it "C+" because it was
originally written in Fortran, then converted to C, and now it's somewhere
between C and C++.  You'll need a C++ compiler.  I tried to structure it as
similarly as I could to the Matlab 2.00 version.  Because it uses a different
method to compute the pseudoinverse, it produces results that are comparable
to, but not identical to, the Matlab version.

Any volunteers for beta testing?  Contact me by removing the "i" from
"comicast".

Greg Berchin