Yeah, while I might be able to make a comment or two, I doubt that I
would be able to pull off the "intelligent" part :-)
Tuomas forgot to attach the patch to his first mail, so I'm attaching it
now.
Thanks,
Søren
man, 18 10 2010 kl. 06:11 -0700, skrev Robert T. Short:
> I have some expertise here. Not in this specific area, but I can
> certainly make some intelligent comments.
>
> Bob
>
> Søren Hauberg wrote:
> > Hi
> >
> > I'm sorry about the long reply.
> >
> > man, 13 09 2010 kl. 01:05 +0300, skrev Tuomas Pennanen:
> >
> >> I found out that the current implementation of interp.m suffers from
> >> noticeable intersymbol interference because of the FIR filter used after
> >> upsampling. I patched interp() to use intfilt() instead of fir1() as per
> >> MathWorks documentation, which solves the problem.
> >>
> >> I also created new file intfilt.m, which is needed by the patched
> >> function. The Lagrange polynomial interpolation in intfilt.m is
> >> unimplemented as I don't know how to do it right.
> >>
> > Is anybody on the list qualified to determine if this is the right
> > change? I don't have the domain specific knowledge to say something
> > qualified...
> >
> >
> > Søren
> >
> >
> > ------------------------------------------------------------------------------
> > Download new Adobe(R) Flash(R) Builder(TM) 4
> > The new Adobe(R) Flex(R) 4 and Flash(R) Builder(TM) 4 (formerly
> > Flex(R) Builder(TM)) enable the development of rich applications that run
> > across multiple browsers and platforms. Download your free trials today!
> > http://p.sf.net/sfu/adobe-dev2dev
> > _______________________________________________
> > Octave-dev mailing list
> > [email protected]
> > https://lists.sourceforge.net/lists/listinfo/octave-dev
> >
>
Index: main/signal/inst/interp.m
===================================================================
--- main/signal/inst/interp.m (revision 7703)
+++ main/signal/inst/interp.m (working copy)
@@ -15,8 +15,11 @@
## usage: y = interp(x, q [, n [, Wc]])
##
-## Upsample the signal x by a factor of q, using an order 2*q*n+1 FIR
-## filter. Note that q must be an integer for this rate change method.
+## Upsample the signal x by a factor of q and calculate the interpolated
+## values using an order 2*q*n-1 FIR filter. The interpolated samples
+## are calculated using the filter returned by intfilt(q, n, Wc).
+##
+## Note that q must be an integer for this rate change method.
## n defaults to 4 and Wc defaults to 0.5.
##
## Example
@@ -27,30 +30,44 @@
## stem(t(1:121)*1000,y(1:121),"-r;Interpolated;");
## stem(t(1:4:121)*1000,x(1:4:121),"-b;Subsampled;"); hold off;
##
-## See also: decimate, resample
+## See also: decimate, intfilt, resample
+## Author: Paul Kienzle
+## Modified by: Tuomas Pennanen, <[email protected]> Sep, 2010
+
function y = interp(x, q, n, Wc)
- if nargin < 1 || nargin > 4,
- usage("y=interp(x, q [, n [, Wc]])");
- endif
- if q != fix(q), error("decimate only works with integer q."); endif
+ if nargin < 1 || nargin > 4
+ usage("y = interp(x, q [, n [, Wc]])");
+ end
- if nargin<3
- n=4; Wc=0.5;
- elseif nargin<4
+ if q ~= fix(q)
+ error("decimate only works with integer q.");
+ end
+
+ if nargin < 3
+ n=4;
+ end
+
+ if nargin < 4
Wc=0.5;
- endif
+ end
- if rows(x)>1
- y = zeros(length(x)*q+q*n+1,1);
+ samples = (length(x) * q) + (q*n - 1);
+
+ if rows(x) > 1
+ x_int = zeros(samples, 1);
else
- y = zeros(1,length(x)*q+q*n+1);
- endif
- y(1:q:length(x)*q) = x;
- b = fir1(2*q*n+1, Wc/q);
- y=q*fftfilt(b, y);
- y(1:q*n+1) = []; # adjust for zero filter delay
+ x_int = zeros(1, samples);
+ end
+ x_int(1:q:(length(x)*q)) = x;
+
+ b = intfilt(q, n, Wc);
+
+ y = fftfilt(b, x_int);
+
+ y(1:(q*n-1)) = []; % adjust for zero filter delay
+
endfunction
%!demo
Index: main/signal/inst/intfilt.m
===================================================================
--- main/signal/inst/intfilt.m (revision 0)
+++ main/signal/inst/intfilt.m (revision 0)
@@ -0,0 +1,58 @@
+% Copyright (C) 2010 Tuomas Pennanen <[email protected]>
+%
+% This program is free software; you can redistribute it and/or modify
+% it under the terms of the GNU General Public License as published by
+% the Free Software Foundation; either version 3 of the License, or
+% (at your option) any later version.
+%
+% This program is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+% GNU General Public License for more details.
+%
+% You should have received a copy of the GNU General Public License
+% along with this program; if not, see <http://www.gnu.org/licenses/>.
+
+% INTFILT(L, P, ALPHA) designs an interpolation filter.
+% B = INTFILT(L, P, ALPHA) designs a filter which blocks aliasing
+% products that are created when a signal is upsampled. L is the
+% upsampling factor and P is the number of samples in the original
+% waveform that contribute to the calculation of an interpolated sample.
+% intfilter assumes that the original data is band limited below the
+% normalized frequency ALPHA, where 0 < ALPHA < 1.
+
+% Author: Tuomas Pennanen <[email protected]>
+
+
+function b = intfilt(l, p, alpha)
+
+ % Length of the filter
+ N = 2*l*p - 1;
+
+ % Scaled passband edge frequency after interpolation
+ Fp = alpha / l;
+ % Scaled Nyquist frequency after interpolation
+ Fnyquist = 1.0 / l;
+
+ % Calculate band edge frequency pairs
+ % Loop as long as F1 < 1
+ last = floor((1 + Fp) / (2*Fnyquist));
+ F = zeros(1, 2*last + 2);
+ for i = 0 : last
+ F(2*i + 1) = i * (2*Fnyquist) - Fp;
+ F(2*i + 2) = i * (2*Fnyquist) + Fp;
+ end
+ F(1) = 0;
+ if F(end) > 1
+ F(end) = 1;
+ end
+
+ % Fill in the amplitude pairs
+ A = zeros(1, 2*last + 2);
+ % First pair specifies the passband gain
+ A(1) = l;
+ A(2) = l;
+
+ b = firls(N-1, F, A);
+
+endfunction
------------------------------------------------------------------------------
Nokia and AT&T present the 2010 Calling All Innovators-North America contest
Create new apps & games for the Nokia N8 for consumers in U.S. and Canada
$10 million total in prizes - $4M cash, 500 devices, nearly $6M in marketing
Develop with Nokia Qt SDK, Web Runtime, or Java and Publish to Ovi Store
http://p.sf.net/sfu/nokia-dev2dev
_______________________________________________
Octave-dev mailing list
[email protected]
https://lists.sourceforge.net/lists/listinfo/octave-dev
