Upsample/FIR/downsample

The example in the upfirdn() - link below - is as follows.

L = 147; M = 160; % Interpolation/decimation factors.
N = 24*M;
h = fir1(N,1/M,kaiser(N+1,7.8562));
h = L*h; % Passband gain = L
Fs = 48e3; % Original sampling frequency-48kHz
n = 0:10239; % 10240 samples, 0.213 seconds long
x = sin(2*pi*1e3/Fs*n); % Original signal, sinusoid @ 1kHz
y = upfirdn(x,h,L,M); % 9408 samples, still .213 seconds

% Overlay original (48kHz) with resampled
% signal (44.1kHz) in red.

stem(n(1:49)/Fs,x(1:49)); hold on
stem(n(1:45)/(Fs*L/M),y(13:57),'r','filled');
xlabel('Time (sec)');ylabel('Signal value');


I think I understand the factors 147 and 160, but cannot figure out
why N = 24*M?

My understanding of

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

is an FIR of order N, cut-off frequency pi/M, and the given Kaiser
window.  Am I correct?

Thanks in advance.


http://www.mathworks.com/access/helpdesk/help/toolbox/signal/index.html?/access/helpdesk/help/toolbox/signal/upfirdn.html&http://www.google.com/search?hl=en&q=upfirdn

On Fri, 03 Aug 2007 19:38:18 -0700, mpowell_jr@hotmail.com wrote:

>
>
>The example in the upfirdn() - link below - is as follows.
>
>L = 147; M = 160; % Interpolation/decimation factors.
>N = 24*M;
>h = fir1(N,1/M,kaiser(N+1,7.8562));
>h = L*h; % Passband gain = L
>Fs = 48e3; % Original sampling frequency-48kHz
>n = 0:10239; % 10240 samples, 0.213 seconds long
>x = sin(2*pi*1e3/Fs*n); % Original signal, sinusoid @ 1kHz
>y = upfirdn(x,h,L,M); % 9408 samples, still .213 seconds
>
>% Overlay original (48kHz) with resampled
>% signal (44.1kHz) in red.
>
>stem(n(1:49)/Fs,x(1:49)); hold on
>stem(n(1:45)/(Fs*L/M),y(13:57),'r','filled');
>xlabel('Time (sec)');ylabel('Signal value');
>
>
>I think I understand the factors 147 and 160, but cannot figure out
>why N = 24*M?
>
>My understanding of
>
>h = fir1(N,1/M,kaiser(N+1,7.8562))
>
>is an FIR of order N, cut-off frequency pi/M, and the given Kaiser
>window.  Am I correct?
>
>Thanks in advance.
>http://www.mathworks.com/access/helpdesk/help/toolbox/signal/index.html?/access/helpdesk/help/toolbox/signal/upfirdn.html&http://www.google.com/search?hl=en&q=upfirdn

Hello mpowell,

Here are my two cents:

Because your interpolation factor is so 
darned high, the filter must be a super 
narrowband (with a super-narrow transition 
region width) filter.   Your desired filter 
passband width is roughly Fs/(2*160)!!
And I'm guessing your filter's transition 
region width is on the same order of magnitude.
Such a filter, having 50-60 dB of stopband 
attenuation, will have many hundreds, or maybe 
a couple thousand taps.  That's why I think 
N was defined to be so large (N = 3840).

Now MATLAB's "upfirdn" command (written by the 
great James McClellan of Parks-McClellan filter 
design fame) is a MEX file so I couldn't learn 
exactly how the filtering is implemented.  But 
the filter is implemented using polyphase 
filter decomposition, so that tells me that 
the filter's tap length mostly likely needs to 
be an integer multiple of L or M (I don't 
know which).  So I'm guessing the polyphase 
implementation is why the filter length 
was set to an integer multiple of M.

BUT!!  Here's the part that puzzles me.
When we use:

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

the number of coefficients in h is
N+1!!  That means the number of taps, 3841, 
is not an integer multiple of either L = 147
or M = 160.  I sure can't explain that.

Shoot mpowell, I don't think I was of much 
help here.  Let's hope somebody else can be.

See Ya',
[-Rick-]

>
> Shoot mpowell, I don't think I was of much
> help here.  Let's hope somebody else can be.
>
> See Ya',
> [-Rick-]- Hide quoted text -

Rick as always I appreaciate the response.  Still soaking up your
comments. If only I knew half as much as you, Jerry, Bristow among
others

mpowell_jr@hotmail.com wrote:
>> Shoot mpowell, I don't think I was of much
>> help here.  Let's hope somebody else can be.
>>
>> See Ya',
>> [-Rick-]- Hide quoted text -
> 
> Rick as always I appreaciate the response.  Still soaking up your
> comments. If only I knew half as much as you, Jerry, Bristow among
> others

Who, me? I don't know nuttin!

Jerry
-- 
Engineering is the art of making what you want from things you can get.
¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

On Aug 3, 10:38 pm, mpowell...@hotmail.com wrote:
> The example in the upfirdn() - link below - is as follows.
>
> L = 147; M = 160; % Interpolation/decimation factors.
> N = 24*M;
> h = fir1(N,1/M,kaiser(N+1,7.8562));
> h = L*h; % Passband gain = L
> Fs = 48e3; % Original sampling frequency-48kHz
> n = 0:10239; % 10240 samples, 0.213 seconds long
> x = sin(2*pi*1e3/Fs*n); % Original signal, sinusoid @ 1kHz
> y = upfirdn(x,h,L,M); % 9408 samples, still .213 seconds
>
> % Overlay original (48kHz) with resampled
> % signal (44.1kHz) in red.
>
> stem(n(1:49)/Fs,x(1:49)); hold on
> stem(n(1:45)/(Fs*L/M),y(13:57),'r','filled');
> xlabel('Time (sec)');ylabel('Signal value');
>
> 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.

essentially, the way they are presenting this problem called:

  "Change the sampling rate by a factor of 147/160. This factor is
used to
  convert from 48kHz (DAT rate) to 44.1kHz (CD sampling rate)."

is that you are conceptually upsampling by a factor of 147 (so 48 kHz
* 147 = 7.056 Mhz) and then downsample by a factor of 160 (7.058 Mhz =
44.1 kHz * 160).  you're doing this (again *conceptually*, you're not
actually gonna bloody your fists doing exactly this - you'll win this
battle with a little subterfuge) by padding 146 zeros between each of
your actual samples.  then, to downsample, you're gonna throw away 159
adjacent samples (that you've fought so hard to compute) out of every
160 and what you got left (after lotsa blood and bruises) is your 44.1
kHz audio.

that's upsampled to 7 Mhz, *but* you still have images at every
multiple of 48 kHz.  one at 48, another at 96, another at -48, etc.
if you were to resample this at 44.1 kHz, all of those images would
fold over into your baseband and you would have an icky mess
resembling canine fecal matter.

so now (**conceptually**, not really) at the sampling rate of 7.056
Mhz, you're gonna low-pass filter them images (except for the original
image centered at 0 Hz) out of the picture with a big, expensive
24*160 = 3840 tap FIR brickwall LPF that gets designed by fir1(.).
the cutoff frequency is gonna be Nyquist (3.528 MHz) divided by 160 or
22.05 kHz.  hence the 1/M in fir1().

but since 146 out of every 147 samples are zero, your big, expensive 7
Mhz FIR is wasting a lot of multiply-accumulate instructions on terms
that you know, *in*advance* are friggin zero.  why bother multiplying
a bunch of pairs of numbers where one of each pair are known to be
zero?  i now see a problem in MATLAB's example.  they should have made
their FIR to be 24*147 = 3528 taps long.  then you would know that out
of the 3528 taps, exactly 24 equally spaced taps will have non-zero
data and the other 3504 taps will have the zero samples.  so their
first mistake is that

h = fir1(24*160,1/M,kaiser(N+1,7.8562));

should be

h = fir1(24*147,1/M,kaiser(N+1,7.8562));

so let's correct that error and proceed as if M = 24*147 = 3528
instead of 3840.  (i'll bet the arrogant assholes at MathWorks do not
admit their error.  they never admit their errors, until MATLAB can
have non-positive indices in their arrays.)

so we know that exactly 24 equally spaced taps have non-zero data and
the opther 3504 have zero, so, to save a little blood, we won't bother
crunching those numbers.  now there are 147 different possible
alignments with the upsampled signal relative to the non-zero
samples.  we are incrementing our index by 160 for each output
sample.  so let's say we start out at 0 and the taps with non-zero
data are at 0, 147, 294, etc.  those are the coefficients we'll use
and the other coefficients we'll skip over.

then we add 160 to our index for our next output sample (we save more
of our precious blood by choosing not to fight the other 159 battles
between tap 0 and tap 160 - fuck 'em, we don't need those FIR outputs
in our final 44.1 kHz signal, we'd be throwing them away anyway).
that scoots us ahead by 160 samples (at the 7 Mhz rate).  the taps at
0 and 147 are behind us but we still have 294, 441, 588... and,
relative to the output sample location of 160, those taps now have
relative indices 134, 281, 428... in the big 3528 tap FIR.

now if you look at this from the POV of the original samples (at 48
kHz), all of the samples are non-zero but our output pointer or index
is not an integer, but also has a fractional value.  for every output
sample, this index advances by a step of 160/147 (which is
48000/44100).  the integer part of this index will tell you which 24
samples to used (all of them at integer indices of equal or ahead of
that integer part).  the fractional part which is one of 147 possible
values, (0/147, 1/147, 2/147 ... 146/147) tells you which set of 24
coefficients out of the 3528 coefs, that you actually have to use to
combine in an FIR filter to get your output for one sample at 44.1
kHz.

if this is not clear, there is something at DSPguru that explains it
this way, but i also periodically post a little FAQ on the sampling
theorem and how to think of resampling from that POV (rather than this
"upsample by 147 and downsample by 160" thing which i consider to be
less intuitive).  i like to think of resampling as hypothetically
reconstructing the continuous-time output (using what we know from the
Nyquist/Shannon Sampling and Reconstruction Theorem) and then sampling
that continuous-time function at the *new* sampling instances.  if
your reconstruction filter (with a continuous-time impulse response)
has a finite length (rather than the perfect sinc() function), it
eventually comes out to be the same math as the above way of looking
at it.

there are a lot of EE topics where i choose not to look at it the same
way as most other EEs, and this is one of them.

> My understanding of
>
> h = fir1(N,1/M,kaiser(N+1,7.8562))
>
> is an FIR of order N, cut-off frequency pi/M, and the given Kaiser
> window.  Am I correct?

yes, almost.

the FIR is of order N-1, but the number of taps is the FIR order + 1
(or N).  but they fucked it up.  N should be 24*147, not 24*160.

r b-j

On Aug 5, 1:02 am, robert bristow-johnson <r...@audioimagination.com>
wrote:
...
>  i now see a problem in MATLAB's example.  they should have made
> their FIR to be 24*147 = 3528 taps long.  then you would know that out
> of the 3528 taps, exactly 24 equally spaced taps will have non-zero
> data and the other 3504 taps will have the zero samples.  so their
> first mistake is that
>
> h = fir1(24*160,1/M,kaiser(N+1,7.8562));
>
> should be
>
> h = fir1(24*147,1/M,kaiser(N+1,7.8562));
>

my correction to MATLAB's code needs to be corrected.  they made two
errors.  i'm surprized they could make this work at all, because
MATLAB is (understandably) quite unforgiving when dimensions of rows
and columns of multiplied matrices are incompatible.  the program
should look like:

L = 147; M = 160;     % Interpolation/decimation factors.
N = 24*L;   % 24 taps out of N will be non-zero for every output
sample.
h = fir1(N-1,1/M,kaiser(N,7.8562));  % FIR order is N-1, but N is the
number of taps
h = L*h; % Passband gain = L
Fs = 48e3;            % Original sampling frequency-48kHz
n = 0:10239;          % 10240 samples, 0.213 seconds long
x  = sin(2*pi*1e3/Fs*n); % Original signal, sinusoid @ 1kHz
y = upfirdn(x,h,L,M);    % 9408 samples, still .213 seconds

% Overlay original (48kHz) with resampled
% signal (44.1kHz) in red.

stem(n(1:49)/Fs,x(1:49)); hold on
stem(n(1:45)/(Fs*L/M),y(13:57),'r','filled');
xlabel('Time (sec)');ylabel('Signal value');

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

r b-j

"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

-- 
Steve Lord
slord@mathworks.com 

On Aug 6, 1:36 pm, "Steven Lord" <sl...@mathworks.com> wrote:
> "robert bristow-johnson" <r...@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/

i needed to register or something that made it painful to report.

> 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/27b53982eae...

give my best to Cleve.  please tell him that it's still not to late to
fix the horrible indexing problem for MATLAB (when using an FFT, DC is
not in "bin #1").

r b-j

On Sun, 05 Aug 2007 05:22:18 -0000, robert bristow-johnson
<rbj@audioimagination.com> wrote:

>On Aug 5, 1:02 am, robert bristow-johnson <r...@audioimagination.com>
>wrote:
>...
>>  i now see a problem in MATLAB's example.  they should have made
>> their FIR to be 24*147 = 3528 taps long.  then you would know that out
>> of the 3528 taps, exactly 24 equally spaced taps will have non-zero
>> data and the other 3504 taps will have the zero samples.  so their
>> first mistake is that
>>
>> h = fir1(24*160,1/M,kaiser(N+1,7.8562));
>>
>> should be
>>
>> h = fir1(24*147,1/M,kaiser(N+1,7.8562));
>>
>
>my correction to MATLAB's code needs to be corrected.  they made two
>errors.  i'm surprized they could make this work at all, because
>MATLAB is (understandably) quite unforgiving when dimensions of rows
>and columns of multiplied matrices are incompatible.  the program
>should look like:
>
>L = 147; M = 160;     % Interpolation/decimation factors.
>N = 24*L;   % 24 taps out of N will be non-zero for every output
>sample.
>h = fir1(N-1,1/M,kaiser(N,7.8562));  % FIR order is N-1, but N is the
>number of taps
>h = L*h; % Passband gain = L
>Fs = 48e3;            % Original sampling frequency-48kHz
>n = 0:10239;          % 10240 samples, 0.213 seconds long
>x  = sin(2*pi*1e3/Fs*n); % Original signal, sinusoid @ 1kHz
>y = upfirdn(x,h,L,M);    % 9408 samples, still .213 seconds
>
>% Overlay original (48kHz) with resampled
>% signal (44.1kHz) in red.
>
>stem(n(1:49)/Fs,x(1:49)); hold on
>stem(n(1:45)/(Fs*L/M),y(13:57),'r','filled');
>xlabel('Time (sec)');ylabel('Signal value');

Hi Robert,
   I'll bet that you're correct.

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 L=160, or M=147.

See Ya',
[-Rick-]