Upsample/FIR/downsample

On Mon, 6 Aug 2007 13:36:36 -0400, "Steven Lord" <slord@mathworks.com>
wrote:

>
  (snipped)
>
>I've asked the developers responsible for the Signal Processing Toolbox to 
>take a look at that example and correct it and/or add some comments 
>clarifying what it's doing.  Thanks for bringing it to our attention.
>
>If you find something like this in the future, please report it to The 
>MathWorks directly in addition to posting it to comp.dsp.  This will make 
>certain that we're aware of the issue (sometimes Usenet messages get delayed 
>or lost) and will let technical support track how many people have reported 
>the issue.  You can send feedback in using the forms on this page:
>
>http://www.mathworks.com/company/aboutus/contact_us/
>
>For those CSSM readers who want to mosey on over to comp.dsp to read the 
>first messages in this thread, the Google groups link is:
>
>http://groups.google.com/group/comp.dsp/browse_frm/thread/27b53982eaed9291/4c4d7ef1c9f78a68?hl=en#4c4d7ef1c9f78a68
>
>-- 
>Steve Lord
>slord@mathworks.com 

Hello Steve,

   Thanks for your trouble.

[-Rick-]
If you run across Ric Losada there at MathWorks
please say "Hello" to him for me.  Thanks.

On Aug 6, 6:04 pm, R.Lyons@_BOGUS_ieee.org (Rick Lyons) wrote:
>
> By the way, in the Math Works "Signal Processing
> ToolboxUser's Guide" textbook, Version Four
> (dated 1996), on page 6-301,
> their upfirdn() example is interpolation by
> a factor of 160/147.  There they designed the FIR
> filter using:
>
>    h = fir1(300,2/160);
>
> So in that textbook example the FIR filter length
> was unrelated to either M=160, or L=147.

i guess i don understand upfirdn(), but my understand of doing sample
rate conversion by a fixed and rational ratio of L/M is that you first
hypothetically upsample by L, then decimate by factor M, both
integers.  so in your hypothetical L-upsampled signal, you have L-1
zeros padded between each original sample at indices L*n (n are the
original indices of x[n]).  that means L different phases and that
means your LPF (with cutoff at min(Nyquist/L, Nyquist/M) ) is really
some N*L samples long and you have L different sets of N sample
impulse responses.  even if they zero-pad what comes out of fir1(),
they must have L different sets of impulse responses of some fixed N
samples in length .

so maybe upfirdn, pads it to be of that length, before applying
reshape() to the big impulse response.

if i were to do this without upfirdn(), it might look like:

n_taps = 24;
L = 147; M = 160;     % Interpolation/decimation factors.
y = zeros(1, floor(length(x)*L/M));

N = n_taps*L;
h = fir1(N-1, min(1/M,1/L), kaiser(N,7.8562));
h = L*h; % Passband gain = L
h = reshape(h, [L n_taps]);

x = [x zeros(1, n_taps)];  %pad a little
for k = 0 : length(y)-1,


   n = floor((M*k)/L);
   phase = M*k - L*n;  % 0 <= phase < L
   y(k+1) = h(phase+1,:) * x(n+1:n+n_taps);
           % dot product
end;

i  dunno  if  i lined everything up right.  i have Octave here at
work, but not MATLAB and am not sure that Octave has stuff like
fir1(.) to design filters.   so this is not tested.  note that my
indices start at zero in the for loop, but i had to add 1 when using
them as a MATLAB index.  now ask what the code would look like if my
for loop started at 1?  it always makes better sense to index from 0.
at least when doing DSP.  it's an irony that probably DSP is the most
widely used application of MATLAB, but they force us to accept that DC
is in bin #1 of the FFT output.


r b-j

Steven Lord wrote:
> "robert bristow-johnson" <rbj@audioimagination.com> wrote in message 
> news:1186291696.843200.210150@w3g2000hsg.googlegroups.com...
>> On Aug 5, 1:02 am, robert bristow-johnson <r...@audioimagination.com>
>> wrote:
>>> On Aug 3, 10:38 pm, mpowell...@hotmail.com wrote:
>>>
>> ...
>>>> I think I understand the factors 147 and 160, but cannot figure out
>>>> why N = 24*M?
>>> i believe that N=24, means that your polyphase filter (for resampling)
>>> is a 24-tap FIR filter and there are 160 different "phases" or
>>> fractional-sample delays for it.
>> and i meant to correct my other misleading statement (following
>> MathWorks' mis-lead), to say that there are 147 different phases or
>> fractional-sample delays, not 160.  i gotta be more careful, even if
>> TMW is not.
>>
>> (now cross-posted to comp.soft-sys.matlab.  you c.s-s.m folks will
>> have to mosey on over to comp.dsp to see what the rest of this is
>> about.  lotsa fun when TMW publishes mistaken code examples.)
> 
> I've asked the developers responsible for the Signal Processing Toolbox to 
> take a look at that example and correct it and/or add some comments 
> clarifying what it's doing.  Thanks for bringing it to our attention.
> 
> If you find something like this in the future, please report it to The 
> MathWorks directly in addition to posting it to comp.dsp.  This will make 
> certain that we're aware of the issue (sometimes Usenet messages get delayed 
> or lost) and will let technical support track how many people have reported 
> the issue.  You can send feedback in using the forms on this page:
> 
> http://www.mathworks.com/company/aboutus/contact_us/
> 
> For those CSSM readers who want to mosey on over to comp.dsp to read the 
> first messages in this thread, the Google groups link is:
> 
> http://groups.google.com/group/comp.dsp/browse_frm/thread/27b53982eaed9291/4c4d7ef1c9f78a68?hl=en#4c4d7ef1c9f78a68
> 

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

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. Using a windowed sinc function eliminates these problems.

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.

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.

Back to the original question - the 24*M factor is just a heuristic
way of estimating the order of the necessary filter.

So far as I know, the kaiser window design is the only window in which
you can estimate the order of the filter for based on requirements.
With the other windows it is more just a case of trial and error.

r b-j  - next you'll be asking for a function that returns the length
of a vector minus 1 :)

Cheers,
David

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 ...

> 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.

> 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.

> 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.

> So far as I know, the kaiser window design is the only window in which
> you can estimate the order of the filter for based on requirements.

no, they got formulae to estimate order for P-McC (in MATLAB, they
used to be called remezord()).

> With the other windows it is more just a case of trial and error.

oh, i see what you mean.  yeah, i think that Kaiser is the only window
that let's you define the stopband attenuation (but not the passband
ripple) and the transistion band and gives you an estimate of the
length.  this is in O&S, i think, but also in the MATLAB description
of kaiser().

> r b-j  - next you'll be asking for a function that returns the length
> of a vector minus 1 :)

:-)

r b-j

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

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, 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 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 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