[FFmpeg-devel] [RFC] Lowpass filter
Michael Niedermayer
michaelni
Tue Aug 26 14:18:18 CEST 2008
On Tue, Aug 26, 2008 at 08:48:54AM +0300, Kostya wrote:
> On Mon, Aug 25, 2008 at 07:13:35PM +0200, Michael Niedermayer wrote:
> > On Mon, Aug 25, 2008 at 12:38:31PM +0300, Kostya wrote:
> > > Here's my lowpass filter made more generic, so it will create
> > > lowpass filter for any even order with any cutoff frequency.
> > >
> > > I submit this in a plain form since it's easier to review
> > > this way (there's not too much left from the old implementation).
> >
> > [...]
> > > #include "avcodec.h"
> > >
> > > struct FFLPFilterCoeffs;
> > > struct FFLPFilterState;
> > >
> > > /**
> > > * Initialize filter coefficients.
> > > *
> > > * @param order filter order
> > > * @param cutoff_ratio cutoff to input frequency ratio
> > > *
> > > * @return pointer to filter coefficients structure or NULL if filter cannot be created
> > > */
> > > struct FFLPFilterCoeffs* ff_lowpass_filter_init_coeffs(int order, float cutoff_ratio);
> > >
> > > /**
> > > * Create new filter state.
> > > *
> > > * @param order filter order
> > > *
> > > * @return pointer to new filter state or NULL if state creation fails
> > > */
> > > struct FFLPFilterState* ff_lowpass_filter_init_state(int order);
> > >
> > > /**
> > > * Free filter coefficients.
> > > *
> > > * @param coeffs pointer allocated with ff_lowpass_filter_init_coeffs()
> > > */
> > > void ff_lowpass_filter_free_coeffs(struct FFLPFilterCoeffs *coeffs);
> > >
> > > /**
> > > * Free filter state.
> > > *
> > > * @param state pointer allocated with ff_lowpass_filter_init_state()
> > > */
> > > void ff_lowpass_filter_free_state(struct FFLPFilterState *state);
> > >
> > > /**
> > > * Perform lowpass filtering on input samples.
> > > *
> > > * @param coeffs pointer to filter coefficients
> > > * @param state pointer to filter state
> > > * @param size input length
> > > * @param src source samples
> > > * @param sstep source stride
> > > * @param dst filtered samples (destination may be the same as input)
> > > * @param dstep destination stride
> > > */
> > > void ff_lowpass_filter(const struct FFLPFilterCoeffs *coeffs, struct FFLPFilterState *state, int size, int16_t *src, int sstep, int16_t *dst, int dstep);
> > >
> > > #endif /* FFMPEG_LOWPASS_H */
> >
> > [...]
> >
> > > struct FFLPFilterCoeffs* ff_lowpass_filter_init_coeffs(int order, float cutoff_ratio)
> > > {
> > > int i, j, size;
> > > FFLPFilterCoeffs *c;
> > > double wa;
> > > complex *p, *zp;
> > >
> > > if(order <= 1 || (order & 1) || cutoff_ratio >= 1.0)
> > > return NULL;
> > >
> > > c = av_malloc(sizeof(FFLPFilterCoeffs));
> > > c->cx = av_mallocz(sizeof(c->cx[0]) * (order + 1));
> > > c->cy = av_malloc (sizeof(c->cy[0]) * order);
> > > c->order = order;
> > >
> > > p = av_malloc(sizeof(p[0]) * (order + 1));
> > > zp = av_malloc(sizeof(zp[0]) * (order + 1));
> >
> > Is there some largest order that makes sense?
> > of so and its not too large the mallocs could be avoided.
>
> Original mkfilter program set maximum order to 10, and I limit
> maximum order to 30 now.
>
> > >
> > > wa = 2 * tan(M_PI * cutoff_ratio / 2.0);
> >
> > M_PI * 0.5 * ...
> > so gcc can multiply the constants out at compile time
> >
> >
> > >
> > > for(i = 0; i <= order; i++){
> > > complex t;
> > > double th = (i + order/2 + 0.5) * M_PI / order;
> > >
> >
> > > t = (cos(th) + I*sin(th)) * wa;
> >
> > t= cexp(I*th)*wa
> >
> >
> > > zp[i] = (2.0 + t) / (2.0 - t);
> > > }
> > >
> >
> > > c->cx[0] = 1;
> >
> > > for(i = 0; i <= order; i++)
> > > for(j = i; j >= 1; j--)
> > > c->cx[j] += c->cx[j - 1];
> >
> > for(i = 1; i <= order; i++)
> > c->cx[i]= c->cx[i-1]*(order-i+1)/i;
> >
> >
> >
> > >
> > > p[0] = 1.0;
> > > for(i = 1; i <= order; i++)
> > > p[i] = 0.0;
> > > for(i = 0; i < order; i++){
> > > for(j = order; j >= 1; j--)
> > > p[j] = -zp[i]*p[j] + p[j - 1];
> > > p[0] *= -zp[i];
> > > }
> >
> > this can be merged into the zp generating loop above and the zp array can
> > be droped
> >
> >
> > > c->gain = creal(p[order]);
> > > for(i = 0; i < order; i++){
> >
> > > complex t = -p[i] / p[order];
> > > c->gain += creal(p[i]);
> > > c->cy[i] = creal(t);
> >
> > t seems like a redundant local varible
> >
> >
> > [...]
> > --
> > Michael GnuPG fingerprint: 9FF2128B147EF6730BADF133611EC787040B0FAB
> >
> > The misfortune of the wise is better than the prosperity of the fool.
> > -- Epicurus
> /*
> * Lowpass IIR filter
Is there anything specific to lowpassing in the code?
If not than IIR filter would be a better name.
[...]
> /**
> * Initialize filter coefficients.
> *
> * @param order filter order
> * @param cutoff_ratio cutoff to input frequency ratio
> *
> * @return pointer to filter coefficients structure or NULL if filter cannot be created
> */
> struct FFLPFilterCoeffs* ff_lowpass_filter_init_coeffs(int order, float cutoff_ratio);
This should also have a enum for filter type.
[...]
> //maximum supported filter order
> #define MAXORDER 30
does it work with order 30 or are roundoff errors causing problems already,
i mean at some order they will ...
>
> struct FFLPFilterCoeffs* ff_lowpass_filter_init_coeffs(int order, float cutoff_ratio)
> {
> int i, j, size;
> FFLPFilterCoeffs *c;
> double wa;
> complex p[MAXORDER + 1];
>
> if(order <= 1 || (order & 1) || order > MAXORDER || cutoff_ratio >= 1.0)
> return NULL;
>
> c = av_malloc(sizeof(FFLPFilterCoeffs));
> c->cx = av_mallocz(sizeof(c->cx[0]) * (order / 2 + 1));
> c->cy = av_malloc (sizeof(c->cy[0]) * order);
> c->order = order;
>
> wa = 2 * tan(M_PI * 0.5 * cutoff_ratio);
>
> c->cx[0] = 1;
> for(i = 1; i < order / 2 + 1; i++)
> c->cx[i] = c->cx[i - 1] * (order - i + 1LL) / i;
>
> p[0] = 1.0;
> for(i = 1; i <= order; i++)
> p[i] = 0.0;
> for(i = 0; i < order; i++){
> complex zp;
> double th = (i + order/2 + 0.5) * M_PI / order;
> zp = cexp(I*th) * wa;
> zp = -(2.0 + zp) / (2.0 - zp);
zp = (zp + 2.0) / (zp - 2.0);
>
> for(j = order; j >= 1; j--)
> p[j] = zp*p[j] + p[j - 1];
> p[0] *= zp;
> }
> c->gain = creal(p[order]);
> for(i = 0; i < order; i++){
> c->gain += creal(p[i]);
> c->cy[i] = creal(-p[i] / p[order]);
> }
> c->gain /= 1 << order;
>
> return c;
> }
>
> struct FFLPFilterState* ff_lowpass_filter_init_state(int order)
> {
> FFLPFilterState *s = av_mallocz(sizeof(FFLPFilterState));
> s->x = av_mallocz(sizeof(s->x[0]) * order);
FFLPFilterState *s = av_mallocz(sizeof(FFLPFilterState) + sizeof(s->x[0]) * order);
with appropriate changes to the struct
[...]
--
Michael GnuPG fingerprint: 9FF2128B147EF6730BADF133611EC787040B0FAB
I have never wished to cater to the crowd; for what I know they do not
approve, and what they approve I do not know. -- Epicurus
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 189 bytes
Desc: Digital signature
URL: <http://lists.mplayerhq.hu/pipermail/ffmpeg-devel/attachments/20080826/cc3b4778/attachment.pgp>
More information about the ffmpeg-devel
mailing list