Commit a169f498 authored by Kostya Shishkov's avatar Kostya Shishkov

Add generic IIR filter interface with Butterworth lowpass filter implementation

and remove obsoleted old lowpass filter.

Originally committed as revision 15005 to svn://svn.ffmpeg.org/ffmpeg/trunk
parent 29ca668f
/*
* IIR filter
* Copyright (c) 2008 Konstantin Shishkov
*
* This file is part of FFmpeg.
*
* FFmpeg is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* FFmpeg 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
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with FFmpeg; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
/**
* @file iirfilter.c
* different IIR filters implementation
*/
#include "iirfilter.h"
#include <complex.h>
#include <math.h>
/**
* IIR filter global parameters
*/
typedef struct FFIIRFilterCoeffs{
int order;
float gain;
int *cx;
float *cy;
}FFIIRFilterCoeffs;
/**
* IIR filter state
*/
typedef struct FFIIRFilterState{
float x[1];
}FFIIRFilterState;
/// maximum supported filter order
#define MAXORDER 30
struct FFIIRFilterCoeffs* ff_iir_filter_init_coeffs(enum IIRFilterType filt_type,
enum IIRFilterMode filt_mode,
int order, float cutoff_ratio,
float stopband, float ripple)
{
int i, j, size;
FFIIRFilterCoeffs *c;
double wa;
complex p[MAXORDER + 1];
if(filt_type != FF_FILTER_TYPE_BUTTERWORTH || filt_mode != FF_FILTER_MODE_LOWPASS)
return NULL;
if(order <= 1 || (order & 1) || order > MAXORDER || cutoff_ratio >= 1.0)
return NULL;
c = av_malloc(sizeof(FFIIRFilterCoeffs));
c->cx = av_malloc(sizeof(c->cx[0]) * ((order >> 1) + 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 >> 1) + 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 >> 1) + 0.5) * M_PI / order;
zp = cexp(I*th) * wa;
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 FFIIRFilterState* ff_iir_filter_init_state(int order)
{
FFIIRFilterState* s = av_mallocz(sizeof(FFIIRFilterState) + sizeof(s->x[0]) * (order - 1));
return s;
}
#define FILTER(i0, i1, i2, i3) \
in = *src * c->gain \
+ c->cy[0]*s->x[i0] + c->cy[1]*s->x[i1] \
+ c->cy[2]*s->x[i2] + c->cy[3]*s->x[i3]; \
res = (s->x[i0] + in )*1 \
+ (s->x[i1] + s->x[i3])*4 \
+ s->x[i2] *6; \
*dst = av_clip_int16(lrintf(res)); \
s->x[i0] = in; \
src += sstep; \
dst += dstep; \
void ff_iir_filter(const struct FFIIRFilterCoeffs *c, struct FFIIRFilterState *s, int size, const int16_t *src, int sstep, int16_t *dst, int dstep)
{
int i;
if(c->order == 4){
for(i = 0; i < size; i += 4){
float in, res;
FILTER(0, 1, 2, 3);
FILTER(1, 2, 3, 0);
FILTER(2, 3, 0, 1);
FILTER(3, 0, 1, 2);
}
}else{
for(i = 0; i < size; i++){
int j;
float in, res;
in = *src * c->gain;
for(j = 0; j < c->order; j++)
in += c->cy[j] * s->x[j];
res = s->x[0] + in + s->x[c->order >> 1] * c->cx[c->order >> 1];
for(j = 1; j < c->order >> 1; j++)
res += (s->x[j] + s->x[c->order - j]) * c->cx[j];
for(j = 0; j < c->order - 1; j++)
s->x[j] = s->x[j + 1];
*dst = av_clip_int16(lrintf(res));
s->x[c->order - 1] = in;
src += sstep;
dst += sstep;
}
}
}
void ff_iir_filter_free_state(struct FFIIRFilterState *state)
{
av_free(state);
}
void ff_iir_filter_free_coeffs(struct FFIIRFilterCoeffs *coeffs)
{
if(coeffs){
av_free(coeffs->cx);
av_free(coeffs->cy);
}
av_free(coeffs);
}
/*
* Lowpass IIR filter
* IIR filter
* Copyright (c) 2008 Konstantin Shishkov
*
* This file is part of FFmpeg.
......@@ -20,27 +20,48 @@
*/
/**
* @file lowpass.h
* lowpass filter interface
* @file iirfilter.h
* IIR filter interface
*/
#ifndef FFMPEG_LOWPASS_H
#define FFMPEG_LOWPASS_H
#ifndef FFMPEG_IIRFILTER_H
#define FFMPEG_IIRFILTER_H
#include "avcodec.h"
struct FFLPFilterCoeffs;
struct FFLPFilterState;
struct FFIIRFilterCoeffs;
struct FFIIRFilterState;
enum IIRFilterType{
FF_FILTER_TYPE_BESSEL,
FF_FILTER_TYPE_BUTTERWORTH,
FF_FILTER_TYPE_CHEBYSHEV,
FF_FILTER_TYPE_ELLIPTIC,
};
enum IIRFilterMode{
FF_FILTER_MODE_LOWPASS,
FF_FILTER_MODE_HIGHPASS,
FF_FILTER_MODE_BANDPASS,
FF_FILTER_MODE_BANDSTOP,
};
/**
* Initialize filter coefficients.
*
* @param filt_type filter type (e.g. Butterworth)
* @param filt_mode filter mode (e.g. lowpass)
* @param order filter order
* @param cutoff_ratio cutoff to input frequency ratio
* @param stopband stopband to input frequency ratio (used by bandpass and bandstop filter modes)
* @param ripple ripple factor (used only in Chebyshev filters)
*
* @return pointer to filter coefficients structure or NULL if filter cannot be created
*/
const struct FFLPFilterCoeffs* ff_lowpass_filter_init_coeffs(int order, float cutoff_ratio);
struct FFIIRFilterCoeffs* ff_iir_filter_init_coeffs(enum IIRFilterType filt_type,
enum IIRFilterMode filt_mode,
int order, float cutoff_ratio,
float stopband, float ripple);
/**
* Create new filter state.
......@@ -49,23 +70,21 @@ const struct FFLPFilterCoeffs* ff_lowpass_filter_init_coeffs(int order, float cu
*
* @return pointer to new filter state or NULL if state creation fails
*/
struct FFLPFilterState* ff_lowpass_filter_init_state(int order);
struct FFIIRFilterState* ff_iir_filter_init_state(int order);
#if 0 //enable with arbitrary order filter implementation, use av_free() for filter state only for now
/**
* Free filter coefficients.
*
* @param coeffs pointer allocated with ff_lowpass_filter_init_coeffs()
* @param coeffs pointer allocated with ff_iir_filter_init_coeffs()
*/
void ff_lowpass_filter_free_coeffs(struct FFLPFilterCoeffs *coeffs);
void ff_iir_filter_free_coeffs(struct FFIIRFilterCoeffs *coeffs);
/**
* Free filter state.
*
* @param state pointer allocated with ff_lowpass_filter_init_state()
* @param state pointer allocated with ff_iir_filter_init_state()
*/
void ff_lowpass_filter_free_state(struct FFLPFilterState *state);
#endif
void ff_iir_filter_free_state(struct FFIIRFilterState *state);
/**
* Perform lowpass filtering on input samples.
......@@ -78,6 +97,7 @@ void ff_lowpass_filter_free_state(struct FFLPFilterState *state);
* @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);
void ff_iir_filter(const struct FFIIRFilterCoeffs *coeffs, struct FFIIRFilterState *state,
int size, const int16_t *src, int sstep, int16_t *dst, int dstep);
#endif /* FFMPEG_LOWPASS_H */
#endif /* FFMPEG_IIRFILTER_H */
/*
* Lowpass IIR filter
* Copyright (c) 2008 Konstantin Shishkov
*
* This file is part of FFmpeg.
*
* FFmpeg is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* FFmpeg 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
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with FFmpeg; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
/**
* @file lowpass.c
* lowpass filter implementation
*/
#include "lowpass.h"
/**********************
* TODO:
* support filters with order != 4
* calculate coefficients for filter instead of taking approximate ones from the table
*********************/
/** filter order */
#define LOWPASS_FILTER_ORDER 4
/**
* IIR filter global parameters
*/
typedef struct FFLPFilterCoeffs{
float gain;
float c[LOWPASS_FILTER_ORDER];
}FFLPFilterCoeffs;
/**
* filter data for 4th order IIR lowpass Butterworth filter
*/
static const FFLPFilterCoeffs lp_filter_coeffs[] = {
{ 9.398085e-01, { -0.0176648009, 0.0000000000, -0.4860288221, 0.0000000000 } },
{ 6.816645e-01, { -0.4646665999, -2.2127207402, -3.9912017501, -3.2380429984 } },
{ 4.998150e-01, { -0.2498216698, -1.3392807613, -2.7693097862, -2.6386277439 } },
{ 3.103469e-01, { -0.0965076902, -0.5977763360, -1.4972580903, -1.7740085241 } },
{ 2.346995e-01, { -0.0557639007, -0.3623690447, -1.0304538354, -1.3066051440 } },
{ 1.528432e-01, { -0.0261686639, -0.1473794606, -0.6204721225, -0.6514716536 } },
{ 6.917529e-02, { -0.0202414073, 0.0780167640, -0.5277442247, 0.3631641670 } },
{ 6.178391e-02, { -0.0223681543, 0.1069446609, -0.5615167033, 0.4883976841 } },
{ 5.298685e-02, { -0.0261686639, 0.1473794606, -0.6204721225, 0.6514716536 } },
{ 2.229030e-02, { -0.0647354087, 0.4172275190, -1.1412129810, 1.4320761385 } },
{ 1.693903e-02, { -0.0823177861, 0.5192354923, -1.3444768251, 1.6365345642 } },
{ 7.374053e-03, { -0.1481421788, 0.8650973862, -1.9894244796, 2.1544844308 } },
{ 5.541768e-03, { -0.1742301048, 0.9921936565, -2.2090801108, 2.3024482658 } },
};
/** cutoff ratios for lp_filter_data[] */
static const float lp_cutoff_ratios[] = {
0.5000000000, 0.4535147392, 0.4166666667, 0.3628117914,
0.3333333333, 0.2916666667, 0.2267573696, 0.2187500000,
0.2083333333, 0.1587301587, 0.1458333333, 0.1133786848,
0.1041666667,
};
/**
* IIR filter state
*/
typedef struct FFLPFilterState{
float x[LOWPASS_FILTER_ORDER];
}FFLPFilterState;
const struct FFLPFilterCoeffs* ff_lowpass_filter_init_coeffs(int order, float cutoff_ratio)
{
int i, size;
//we can create only order-4 filters with cutoff ratio <= 0.5 for now
if(order != LOWPASS_FILTER_ORDER) return NULL;
size = sizeof(lp_cutoff_ratios) / sizeof(lp_cutoff_ratios[0]);
if(cutoff_ratio > lp_cutoff_ratios[0])
return NULL;
for(i = 0; i < size; i++){
if(cutoff_ratio >= lp_cutoff_ratios[i])
break;
}
if(i == size)
i = size - 1;
return &lp_filter_coeffs[i];
}
struct FFLPFilterState* ff_lowpass_filter_init_state(int order)
{
if(order != LOWPASS_FILTER_ORDER) return NULL;
return av_mallocz(sizeof(FFLPFilterState));
}
#define FILTER(i0, i1, i2, i3) \
in = *src * c->gain \
+ c->c[0]*s->x[i0] + c->c[1]*s->x[i1] \
+ c->c[2]*s->x[i2] + c->c[3]*s->x[i3]; \
res = (s->x[i0] + in )*1 \
+ (s->x[i1] + s->x[i3])*4 \
+ s->x[i2] *6; \
*dst = av_clip_int16(lrintf(res)); \
s->x[i0] = in; \
src += sstep; \
dst += dstep; \
void ff_lowpass_filter(const struct FFLPFilterCoeffs *c, struct FFLPFilterState *s, int size, int16_t *src, int sstep, int16_t *dst, int dstep)
{
int i;
for(i = 0; i < size; i += 4){
float in, res;
FILTER(0, 1, 2, 3);
FILTER(1, 2, 3, 0);
FILTER(2, 3, 0, 1);
FILTER(3, 0, 1, 2);
}
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment