From 686025404486c751c808aa7441da19cacef637b4 Mon Sep 17 00:00:00 2001 From: Alex Converse Date: Fri, 30 Jan 2009 20:15:48 +0000 Subject: [PATCH] Add the rdft family of transforms (fft/ifft of an all real sequence) to dsputil. Originally committed as revision 16864 to svn://svn.ffmpeg.org/ffmpeg/trunk --- configure | 2 + libavcodec/Makefile | 1 + libavcodec/dsputil.h | 31 +++++++++++ libavcodec/fft.c | 2 +- libavcodec/rdft.c | 127 +++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 162 insertions(+), 1 deletion(-) create mode 100644 libavcodec/rdft.c diff --git a/configure b/configure index fc8529f39b..3f3731ba06 100755 --- a/configure +++ b/configure @@ -102,6 +102,7 @@ show_help(){ echo " --disable-fft disable FFT code" echo " --disable-golomb disable Golomb code" echo " --disable-mdct disable MDCT code" + echo " --disable-rdft disable RDFT code" echo " --enable-hardcoded-tables use hardcoded tables instead of runtime generation" echo " --enable-memalign-hack emulate memalign, interferes with memory debuggers" echo " --enable-beos-netserver enable BeOS netserver" @@ -781,6 +782,7 @@ CONFIG_LIST=" nonfree postproc powerpc_perf + rdft shared small static diff --git a/libavcodec/Makefile b/libavcodec/Makefile index 9abedfc607..7ccf694c2e 100644 --- a/libavcodec/Makefile +++ b/libavcodec/Makefile @@ -29,6 +29,7 @@ OBJS-$(CONFIG_ENCODERS) += faandct.o jfdctfst.o jfdctint.o OBJS-$(CONFIG_FFT) += fft.o OBJS-$(CONFIG_GOLOMB) += golomb.o OBJS-$(CONFIG_MDCT) += mdct.o +OBJS-$(CONFIG_RDFT) += rdft.o OBJS-$(CONFIG_OLDSCALER) += imgresample.o # decoders/encoders diff --git a/libavcodec/dsputil.h b/libavcodec/dsputil.h index 165aef7a09..66557138bf 100644 --- a/libavcodec/dsputil.h +++ b/libavcodec/dsputil.h @@ -674,6 +674,8 @@ typedef struct FFTContext { void (*imdct_half)(struct MDCTContext *s, FFTSample *output, const FFTSample *input); } FFTContext; +extern FFTSample* ff_cos_tabs[13]; + /** * Sets up a complex FFT. * @param nbits log2 of the length of the input array @@ -759,6 +761,35 @@ void ff_imdct_half_sse(MDCTContext *s, FFTSample *output, const FFTSample *input void ff_mdct_calc(MDCTContext *s, FFTSample *out, const FFTSample *input); void ff_mdct_end(MDCTContext *s); +/* Real Discrete Fourier Transform */ + +enum RDFTransformType { + RDFT, + IRDFT, + RIDFT, + IRIDFT, +}; + +typedef struct { + int nbits; + int inverse; + int sign_convention; + + /* pre/post rotation tables */ + FFTSample *tcos; + FFTSample *tsin; + FFTContext fft; +} RDFTContext; + +/** + * Sets up a real FFT. + * @param nbits log2 of the length of the input array + * @param trans the type of transform + */ +int ff_rdft_init(RDFTContext *s, int nbits, enum RDFTransformType trans); +void ff_rdft_calc(RDFTContext *s, FFTSample *data); +void ff_rdft_end(RDFTContext *s); + #define WRAPPER8_16(name8, name16)\ static int name16(void /*MpegEncContext*/ *s, uint8_t *dst, uint8_t *src, int stride, int h){\ return name8(s, dst , src , stride, h)\ diff --git a/libavcodec/fft.c b/libavcodec/fft.c index db79fc90e9..47ecdb48de 100644 --- a/libavcodec/fft.c +++ b/libavcodec/fft.c @@ -42,7 +42,7 @@ DECLARE_ALIGNED_16(FFTSample, ff_cos_8192[4096]); DECLARE_ALIGNED_16(FFTSample, ff_cos_16384[8192]); DECLARE_ALIGNED_16(FFTSample, ff_cos_32768[16384]); DECLARE_ALIGNED_16(FFTSample, ff_cos_65536[32768]); -static FFTSample *ff_cos_tabs[] = { +FFTSample *ff_cos_tabs[] = { ff_cos_16, ff_cos_32, ff_cos_64, ff_cos_128, ff_cos_256, ff_cos_512, ff_cos_1024, ff_cos_2048, ff_cos_4096, ff_cos_8192, ff_cos_16384, ff_cos_32768, ff_cos_65536, }; diff --git a/libavcodec/rdft.c b/libavcodec/rdft.c new file mode 100644 index 0000000000..36c6f1bf05 --- /dev/null +++ b/libavcodec/rdft.c @@ -0,0 +1,127 @@ +/* + * (I)RDFT transforms + * Copyright (c) 2009 Alex Converse + * + * 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 + */ +#include +#include "dsputil.h" + +/** + * @file rdft.c + * (Inverse) Real Discrete Fourier Transforms. + */ + +/* sin(2*pi*x/n) for 0<=xnbits = nbits; + s->inverse = trans == IRDFT || trans == IRIDFT; + s->sign_convention = trans == RIDFT || trans == IRIDFT ? 1 : -1; + + if (nbits < 4 || nbits > 16) + return -1; + + if (ff_fft_init(&s->fft, nbits-1, trans == IRDFT || trans == RIDFT) < 0) + return -1; + + s->tcos = ff_cos_tabs[nbits-4]; + s->tsin = ff_sin_tabs[nbits-4]+(trans == RDFT || trans == IRIDFT)*(n>>2); + for (i = 0; i < (n>>2); i++) { + s->tcos[i] = cos(i*theta); + s->tsin[i] = sin(i*theta); + } + return 0; +} + +/** Map one real FFT into two parallel real even and odd FFTs. Then interleave + * the two real FFTs into one complex FFT. Unmangle the results. + * ref: http://www.engineeringproductivitytools.com/stuff/T0001/PT10.HTM + */ +void ff_rdft_calc_c(RDFTContext* s, FFTSample* data) +{ + int i, i1, i2; + FFTComplex ev, od; + const int n = 1 << s->nbits; + const float k1 = 0.5; + const float k2 = 0.5 - s->inverse; + const FFTSample *tcos = s->tcos; + const FFTSample *tsin = s->tsin; + + if (!s->inverse) { + ff_fft_permute(&s->fft, (FFTComplex*)data); + ff_fft_calc(&s->fft, (FFTComplex*)data); + } + /* i=0 is a special case because of packing, the DC term is real, so we + are going to throw the N/2 term (also real) in with it. */ + ev.re = data[0]; + data[0] = ev.re+data[1]; + data[1] = ev.re-data[1]; + for (i = 1; i < (n>>2); i++) { + i1 = 2*i; + i2 = n-i1; + /* Separate even and odd FFTs */ + ev.re = k1*(data[i1 ]+data[i2 ]); + od.im = -k2*(data[i1 ]-data[i2 ]); + ev.im = k1*(data[i1+1]-data[i2+1]); + od.re = k2*(data[i1+1]+data[i2+1]); + /* Apply twiddle factors to the odd FFT and add to the even FFT */ + data[i1 ] = ev.re + od.re*tcos[i] - od.im*tsin[i]; + data[i1+1] = ev.im + od.im*tcos[i] + od.re*tsin[i]; + data[i2 ] = ev.re - od.re*tcos[i] + od.im*tsin[i]; + data[i2+1] = -ev.im + od.im*tcos[i] + od.re*tsin[i]; + } + data[2*i+1]=s->sign_convention*data[2*i+1]; + if (s->inverse) { + data[0] *= k1; + data[1] *= k1; + ff_fft_permute(&s->fft, (FFTComplex*)data); + ff_fft_calc(&s->fft, (FFTComplex*)data); + } +} + +void ff_rdft_calc(RDFTContext *s, FFTSample *data) +{ + ff_rdft_calc_c(s, data); +} + +av_cold void ff_rdft_end(RDFTContext *s) +{ + ff_fft_end(&s->fft); +}