On Aug 9, 12:22 pm, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote:
> On Wed, 08 Aug 2007 21:28:07 -0000, robert bristow-johnson
>
>
>
> <r...@audioimagination.com> wrote:
> >On Aug 8, 4:21 pm, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote:
> >> On Wed, 08 Aug 2007 17:31:32 -0000, robert bristow-johnson
>
> >> <r...@audioimagination.com> wrote:
>
> >> >it depends on if you are ultimately upsampling (M<L) or downsampling
> >> >(M>L).  that's why i said
>
> >> >    h = fir1(N-1, min(1/M,1/L), kaiser(N,7.8562));
>
> >...
> >>     Regarding the fir1's freq argument, when
> >> I looked at mpowell's original decimation
> >> code, and the interpolation code example in MATLAB's
> >> hardcopy User's Manual, I originally thought
> >> mpowell made an error in his
>
> >>   h = fir1(N,1/M,kaiser(N+1,7.8562))
>
> >> where M = 160.
>
> >> It took me a while to figure it out (and it
> >> surprised me a little), with MATLAB's formating
> >> the fir1's freq argument should always be 1/160
> >> REGARDLESS of whether your interpolating or
> >> decimating!!
>
> >i don't think i agree, Rick.  if you're UPsampling (switch the values
> >for M and L around, M=147 and L=160) then you want the cutoff
> >frequency to be Nyquist/L, no?
>
> >e.g. converting a sound file from CD sampling rate (44.1 kHz) to DAT
> >or DAW or "studio" sampling rate (48 kHz).  so now the intermediate
> >upsampled frequency is the same 7.058 MHz (44.1kHz*160), but now you
> >want your cutoff to be half of the source sampling rate (22.05 kHz),
> >not half of the destination sampling rate (which is higher, 24 kHz).
> >your 7 MHz intermediate signal will have images at 44.1, 88.2 and so
> >on.  there might be something between 22.05 and 24 kHz that you don't
> >like.
>
> >when up-converting, there should be no loss of information, but there
> >will be holes in the spectrum of the intermediate over-sampled
> >signal.  when down-converting, to prevent aliasing, you must filter
> >your input signal to even less of its Nyquist (to that of the final
> >Nyquist).
>
> >i think it should be the Nyquist*min(1/L,1/M), and in fir1() (and all
> >other things MATLAB), "Nyquist"=1.
>
> >that's my story, and i'm stickin to it.
>
> >r b-j
>
> Well, if your text:
>
>    "min(1/L,1/M)"
>
> means "the minimum of 1/147 or 1/160", then
> min(1/L,1/M) always equals 1/160, right?

yah.

> All I was tryin' to say was that the fir1's
> freq argument should always be 1/160 whether
> we're interpolating by 160/147 or decimating
> by 147/160.
>
> As far as I can tell we're in agreement.

yup.  i might have misread or misunderstood what you said earlier.
but the length of the filter should always be N = n_taps*L, so the
design of the filter is a little different in the cases of conversion
ratios, L/M, of 160/147 vs. 147/160.

r b-j

On Wed, 08 Aug 2007 21:28:07 -0000, robert bristow-johnson
<rbj@audioimagination.com> wrote:

>On Aug 8, 4:21 pm, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote:
>> On Wed, 08 Aug 2007 17:31:32 -0000, robert bristow-johnson
>>
>> <r...@audioimagination.com> wrote:
>>
>>
>> >it depends on if you are ultimately upsampling (M<L) or downsampling
>> >(M>L).  that's why i said
>>
>> >    h = fir1(N-1, min(1/M,1/L), kaiser(N,7.8562));
>>
>...
>>     Regarding the fir1's freq argument, when
>> I looked at mpowell's original decimation
>> code, and the interpolation code example in MATLAB's
>> hardcopy User's Manual, I originally thought
>> mpowell made an error in his
>>
>>   h = fir1(N,1/M,kaiser(N+1,7.8562))
>>
>> where M = 160.
>>
>> It took me a while to figure it out (and it
>> surprised me a little), with MATLAB's formating
>> the fir1's freq argument should always be 1/160
>> REGARDLESS of whether your interpolating or
>> decimating!!
>
>i don't think i agree, Rick.  if you're UPsampling (switch the values
>for M and L around, M=147 and L=160) then you want the cutoff
>frequency to be Nyquist/L, no?
>
>e.g. converting a sound file from CD sampling rate (44.1 kHz) to DAT
>or DAW or "studio" sampling rate (48 kHz).  so now the intermediate
>upsampled frequency is the same 7.058 MHz (44.1kHz*160), but now you
>want your cutoff to be half of the source sampling rate (22.05 kHz),
>not half of the destination sampling rate (which is higher, 24 kHz).
>your 7 MHz intermediate signal will have images at 44.1, 88.2 and so
>on.  there might be something between 22.05 and 24 kHz that you don't
>like.
>
>when up-converting, there should be no loss of information, but there
>will be holes in the spectrum of the intermediate over-sampled
>signal.  when down-converting, to prevent aliasing, you must filter
>your input signal to even less of its Nyquist (to that of the final
>Nyquist).
>
>i think it should be the Nyquist*min(1/L,1/M), and in fir1() (and all
>other things MATLAB), "Nyquist"=1.
>
>that's my story, and i'm stickin to it.
>
>r b-j

Well, if your text:

   "min(1/L,1/M)"

means "the minimum of 1/147 or 1/160", then 
min(1/L,1/M) always equals 1/160, right.
All I was tryin' to say was that the fir1's 
freq argument should always be 1/160 whether 
we're interpolating by 160/147 or decimating 
by 147/160.

As far as I can tell we're in agreement.

[-Rick-]

On Aug 9, 2:14 am, "Ron N." <rhnlo...@yahoo.com> wrote:
> On Aug 8, 2:28 pm, robert bristow-johnson <r...@audioimagination.com>
> wrote:
>
> > On Aug 8, 4:21 pm, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote:
>
> > > On Wed, 08 Aug 2007 17:31:32 -0000, robert bristow-johnson
>
> > > <r...@audioimagination.com> wrote:
>
> > > >it depends on if you are ultimately upsampling (M<L) or downsampling
> > > >(M>L).  that's why i said
>
> > > >    h = fir1(N-1, min(1/M,1/L), kaiser(N,7.8562));
>
...
>
> Wouldn't all the reconstruction filters on the
> way between two sample rates have the same cutoff;
> e.g. min(f1/2, f2/2) ?

yes, i think so.  that why i said

     h = fir1(N-1, min(1/M,1/L), kaiser(N,7.8562));

(N = L*n_taps in either case.)

> In the upsampling case one wants to make sure that
> any new samples don't introduce any new image
> spectral content.  And in the downsampling case,
> one wants to low pass filter first so that there
> will be no content at or above Fs/2 of the
> destination sample rate.

i think we're in agreement.

> I prefer to think of resampling as reconstruction,
> with all the polyphase and up/down intermediate rate
> stuff as just (major) implementation optimizations.

in my long treatise on this thread, i tried to put it in those terms.
but i find it easier to think of resampling as emulating, in a purely
digital environment, what happens in the ideal if you were to, at the
first sampling rate, output your samples to an ideal D/A with an ideal
anti-imaging filter, and then resampling at the new sampling instances
spaced apart by the reciprocal of the new sampling rate.  that's the
easier way for me to imagine it and the details of what needs to
happen.

r b-j

On Aug 8, 2:28 pm, robert bristow-johnson <r...@audioimagination.com>
wrote:
> On Aug 8, 4:21 pm, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote:
>
>
>
> > On Wed, 08 Aug 2007 17:31:32 -0000, robert bristow-johnson
>
> > <r...@audioimagination.com> wrote:
>
> > >it depends on if you are ultimately upsampling (M<L) or downsampling
> > >(M>L).  that's why i said
>
> > >    h = fir1(N-1, min(1/M,1/L), kaiser(N,7.8562));
>
> ...
> >     Regarding the fir1's freq argument, when
> > I looked at mpowell's original decimation
> > code, and the interpolation code example in MATLAB's
> > hardcopy User's Manual, I originally thought
> > mpowell made an error in his
>
> >   h = fir1(N,1/M,kaiser(N+1,7.8562))
>
> > where M = 160.
>
> > It took me a while to figure it out (and it
> > surprised me a little), with MATLAB's formating
> > the fir1's freq argument should always be 1/160
> > REGARDLESS of whether your interpolating or
> > decimating!!
>
> i don't think i agree, Rick.  if you're UPsampling (switch the values
> for M and L around, M=147 and L=160) then you want the cutoff
> frequency to be Nyquist/L, no?
>
> e.g. converting a sound file from CD sampling rate (44.1 kHz) to DAT
> or DAW or "studio" sampling rate (48 kHz).

Wouldn't all the reconstruction filters on the
way between two sample rates have the same cutoff;
e.g. min(f1/2, f2/2) ?
In the upsampling case one wants to make sure that
any new samples don't introduce any new image
spectral content.  And in the downsampling case,
one wants to low pass filter first so that there
will be no content at or above Fs/2 of the
destination sample rate.

I prefer to think of resampling as reconstruction,
with all the polyphase and up/down intermediate rate
stuff as just (major) implementation optimizations.


IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M

On Aug 8, 4:21 pm, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote:
> On Wed, 08 Aug 2007 17:31:32 -0000, robert bristow-johnson
>
> <r...@audioimagination.com> wrote:
>
>
> >it depends on if you are ultimately upsampling (M<L) or downsampling
> >(M>L).  that's why i said
>
> >    h = fir1(N-1, min(1/M,1/L), kaiser(N,7.8562));
>
...
>     Regarding the fir1's freq argument, when
> I looked at mpowell's original decimation
> code, and the interpolation code example in MATLAB's
> hardcopy User's Manual, I originally thought
> mpowell made an error in his
>
>   h = fir1(N,1/M,kaiser(N+1,7.8562))
>
> where M = 160.
>
> It took me a while to figure it out (and it
> surprised me a little), with MATLAB's formating
> the fir1's freq argument should always be 1/160
> REGARDLESS of whether your interpolating or
> decimating!!

i don't think i agree, Rick.  if you're UPsampling (switch the values
for M and L around, M=147 and L=160) then you want the cutoff
frequency to be Nyquist/L, no?

e.g. converting a sound file from CD sampling rate (44.1 kHz) to DAT
or DAW or "studio" sampling rate (48 kHz).  so now the intermediate
upsampled frequency is the same 7.058 MHz (44.1kHz*160), but now you
want your cutoff to be half of the source sampling rate (22.05 kHz),
not half of the destination sampling rate (which is higher, 24 kHz).
your 7 MHz intermediate signal will have images at 44.1, 88.2 and so
on.  there might be something between 22.05 and 24 kHz that you don't
like.

when up-converting, there should be no loss of information, but there
will be holes in the spectrum of the intermediate over-sampled
signal.  when down-converting, to prevent aliasing, you must filter
your input signal to even less of its Nyquist (to that of the final
Nyquist).

i think it should be the Nyquist*min(1/L,1/M), and in fir1() (and all
other things MATLAB), "Nyquist"=1.

that's my story, and i'm stickin to it.

r b-j

On Wed, 08 Aug 2007 17:31:32 -0000, robert bristow-johnson
<rbj@audioimagination.com> wrote:

  (snipped)
>
>it depends on if you are ultimately upsampling (M<L) or downsampling
>(M>L).  that's why i said
>
>    h = fir1(N-1, min(1/M,1/L), kaiser(N,7.8562));
>
  (snipped)

Hi,
    Regarding the fir1's freq argument, when 
I looked at mpowell's original decimation 
code, and the interpolation code example in MATLAB's 
hardcopy User's Manual, I originally thought 
mpowell made an error in his 

  h = fir1(N,1/M,kaiser(N+1,7.8562))

where M = 160.

It took me a while to figure it out (and it 
surprised me a little), with MATLAB's formating 
the fir1's freq argument should always be 1/160 
REGARDLESS of whether your interpolating or 
decimating!!

Interesting, huh?
[-Rick-]

On Aug 7, 2:47 pm, dspg...@netscape.net wrote:
> On Aug 7, 12:32 pm, robert bristow-johnson <r...@audioimagination.com>
> wrote:
>
> > On Aug 7, 10:57 am, dspg...@netscape.net wrote:
>
...
> > > Back to the original question - the 24*M factor is just a heuristic
> > > way of estimating the order of the necessary filter.
>
> > 24 is n_taps.  the number of taps in the FIR or dot-product
> > combination that creates one interpolated output sample.  it may be a
> > number that someone pulled out of their butt, but it has meaning.  it
> > could have been 32 (i've seen and done a few different bandlimited
> > interpolators where n_taps is 32), or more.  Bob Adams' original
> > AD1890 SRC chip had n_taps=64.  the number of polyphase FIR
> > taps,n_taps, times the number of polyphases, L, is equal to the length
> > of the big mondo LPF FIR (with cutoff of 1/M) that gets designed by P-
> > McC or firls() or windowed sinc() or whatever.
>
> By heuristic I just meant it was rather arbitrary ... not that it
> didn't have meaning. To me it implies the performance of the overall
> LPF.
>
> If the filter length is determined by the interpolation factor then
> yes, I believe, you are correct it is the length of each polyphase
> filter.  But in this case the filter length should be tied to the
> decimation factor since that is the more stringent requirement,

it depends on if you are ultimately upsampling (M<L) or downsampling
(M>L).  that's why i said

    h = fir1(N-1, min(1/M,1/L), kaiser(N,7.8562));

in my "corrected" version.  if M was smaller than L, then the cutoff
of the brickwall filter should be around Nyquist/L, not Nyquist/M.
that 7.8562 number comes from a formula i remember seeing.  i'll guess
it corresponds to around -80 dB max stopband ripple (or was that
Kaiser factor closer to 5?)

> given
> the numbers in the original code. Vincent's code does the correct
> thing by taking the maximum of the interpolation/decimation values to
> determine the required transition width -

yes, i just didn't know what:

    d = fdesign.rsrc(L,M,'Nyquist',M,TW,Ast);
    H = design(d,'kaiserwin');

 does.  i *do* know what fir1() does.  he's using the max of M or L, i
am doing the equivalent, picking the min of their reciprocals.

> which you would then use
> then use in a remezord type function to determine the required filter
> length. You can still make it fit nicely into a polyphase filter by
> rounding up the filter length to a multiple of the interpolation
> factor to get the most use out the coefficients.

yes, we agree on that.  it makes no sense, IMO, not to do that.

r b-j

On Tue, 07 Aug 2007 08:28:51 -0400, Vincent Pellissier
<vincentp@mathworks.com> wrote:

  (snipped)
>
>Rick and r b-j are correct. upfirdn implements L phases, not M. 
>Therefore it would make more sense to compute the order of the filter as 
>Lp*L-1 where Lp is the length of each polyphase filter. Thanks for 
>reporting this; we will update the example in upfirdn and the associated 
>documentation.
>
>Notice however that the filter in upfirdn example, although not optimal, 
>do the job (i.e. it limits aliasing). It simply uses 27 taps for each 
>polyphase filter. Furthermore, if you have access to Filter Design 
>Toolbox, you can use a minimum order design instead of guessing the 
>filter order. For example:
>
>L = 147; M = 160;     % Interpolation/decimation factors.
>TW = 1/(2*max(L,M));  % Transition Width
>Ast = 80;             % Stopband Attenuation
>d = fdesign.rsrc(L,M,'Nyquist',M,TW,Ast);
>H = design(d,'kaiserwin');
>p = polyphase(H);
>size(p)
>
>ans =
>
>    147    22
>
>A transition width of 0.003125*pi rad/sample and a stopband attenuation 
>of 80 dB yields a polyphase length of 22.
>
>Vincent Pellissier

Hi Vincent,
  Thanks much for the info!

[-Rick-]

On Aug 7, 12:32 pm, robert bristow-johnson <r...@audioimagination.com>
wrote:
> On Aug 7, 10:57 am, dspg...@netscape.net wrote:
>
> > A couple of things to note:
>
> > upfirdn first interpolates with the provided filter and then decimate
> > (without applying any filter). In this case the decimation factor is
> > larger so it actually should determine the filter transition width in
> > order to eliminate aliasing.
>
> > Since the filter is very large, using remez design techniques become
> > numerically unstable. I've found the limit in matlab to be around 2000
> > taps or so. In some cases it will work, but for the majority of times
> > it doesn't.
>
> geez, i couldn't get it to behave for 1024.  if L*n_taps was greater
> than 512, i had to ...

The remezord calculations have difficulties when the band edges are
either near 0 or pi. The handbook formulas don't take into account
where in the frequency band the transition band is  - they only
require the percentage of the overall sampling frequency. You can try
increasing the number of grid points to improve convergence, but this
also slows down the algorithm somewhat. Even given all that there are
still some degenerate cases where the remez algorithm fails to
converge.

>
> > Using a windowed sinc function eliminates these problems.
>
> ... either use a windowed sinc function or (this was messy), used the
> shorter P-McC design to interpolate itself to get a more densely
> sampled impulse response.  kinda sleezy.

Actually you can see this as doing the interpolation in multiple steps
rather than in single step. This effectively allows you to achieve
much higher filter orders. The decimation can also be done in multiple
steps as well.  Unfortunately the implementation / coding tends to get
more complicated.
>
> > The order of the filter is pretty much irrelevant in this case you can
> > use either N or N-1 or N+1. It is really only important when you need
> > particular behaviours at specific frequencies like pi.
>
> no, you want the order of the FIR to be L*n_taps-1 (since the length
> of the FIR is one bigger than the order).  if n_taps and L are both
> integers (and it's hard for me to conceive them not integers), that
> restricts the meaningful choices of the length of the big FIR
> brickwall LPF.
>
> > The polyphase discussion is also irrelevant wrt the original question.
> > Polyphase is more an implementation issue affecting efficiency. The
> > filter can just be padded with zeros to make the length a multiple of
> > the necessary factor.
>
> but, even with zeros padded, your length is still L*n_taps.  it's just
> a crappier LPF (with poorer transition band performance) because the
> effective length is shorter.

Yes, I see what your saying. You can achieve better filtering by
actually using the remaining coefficients rather than zero padding
them.

>
> > Back to the original question - the 24*M factor is just a heuristic
> > way of estimating the order of the necessary filter.
>
> 24 is n_taps.  the number of taps in the FIR or dot-product
> combination that creates one interpolated output sample.  it may be a
> number that someone pulled out of their butt, but it has meaning.  it
> could have been 32 (i've seen and done a few different bandlimited
> interpolators where n_taps is 32), or more.  Bob Adams' original
> AD1890 SRC chip had n_taps=64.  the number of polyphase FIR
> taps,n_taps, times the number of polyphases, L, is equal to the length
> of the big mondo LPF FIR (with cutoff of 1/M) that gets designed by P-
> McC or firls() or windowed sinc() or whatever.

By heuristic I just meant it was rather arbitrary ... not that it
didn't have meaning. To me it implies the performance of the overall
LPF.

If the filter length is determined by the interpolation factor then
yes, I believe, you are correct it is the length of each polyphase
filter.  But in this case the filter length should be tied to the
decimation factor since that is the more stringent requirement, given
the numbers in the original code. Vincent's code does the correct
thing by taking the maximum of the interpolation/decimation values to
determine the required transition width - which you would then use
then use in a remezord type function to determine the required filter
length. You can still make it fit nicely into a polyphase filter by
rounding up the filter length to a multiple of the interpolation
factor to get the most use out the coefficients.

Cheers,
David