ref: 7351e286c81e98b10ddebfbe3e7c31c866fb3ad5
parent: 6211c90def21b649e7c30e75ea1ea07e2a4c8483
author: Jean-Marc Valin <[email protected]>
date: Fri Feb 8 10:09:45 EST 2008
Everything converted to use kiss-fft. Got rid of smallft and fftwrap. Code much leaner now.
--- a/libcelt/Makefile.am
+++ b/libcelt/Makefile.am
@@ -15,9 +15,9 @@
# Sources for compilation in the library
libcelt_la_SOURCES = bands.c bitrdec.c bitree.c bitrenc.c celt.c cwrs.c \
- ecintrin.h entcode.c entdec.c entenc.c fftwrap.c header.c kiss_fft.c \
+ ecintrin.h entcode.c entdec.c entenc.c header.c kiss_fft.c \
kiss_fftr.c laplace.c mdct.c modes.c pitch.c psy.c quant_bands.c \
- quant_pitch.c rangedec.c rangeenc.c rate.c smallft.c vq.c
+ quant_pitch.c rangedec.c rangeenc.c rate.c vq.c
#noinst_HEADERS =
@@ -24,10 +24,10 @@
libcelt_la_LDFLAGS = -version-info @CELT_LT_CURRENT@:@CELT_LT_REVISION@:@CELT_LT_AGE@
noinst_HEADERS = arch.h bands.h bitrdec.h bitree.h bitrenc.h cwrs.h \
- ecintrin.h entcode.h entdec.h entenc.h fftwrap.h kiss_fft.h \
+ ecintrin.h entcode.h entdec.h entenc.h kiss_fft.h \
kiss_fftr.h _kiss_fft_guts.h laplace.h mdct.h \
mfrngcod.h modes.h os_support.h pgain_table.h pitch.h psy.h \
- quant_bands.h quant_pitch.h rate.h smallft.h vq.h
+ quant_bands.h quant_pitch.h rate.h vq.h
noinst_PROGRAMS = testcelt
testcelt_SOURCES = testcelt.c
--- a/libcelt/celt.c
+++ b/libcelt/celt.c
@@ -34,7 +34,7 @@
#include <math.h>
#include "celt.h"
#include "pitch.h"
-#include "fftwrap.h"
+#include "kiss_fftr.h"
#include "bands.h"
#include "modes.h"
#include "entcode.h"
@@ -63,7 +63,7 @@
float *preemph_memD;
mdct_lookup mdct_lookup;
- void *fft;
+ kiss_fftr_cfg fft;
float *window;
float *in_mem;
@@ -97,7 +97,7 @@
ec_enc_init(&st->enc,&st->buf);
mdct_init(&st->mdct_lookup, 2*N);
- st->fft = spx_fft_init(MAX_PERIOD*C);
+ st->fft = kiss_fftr_alloc(MAX_PERIOD*C, 0, 0);
st->window = celt_alloc(2*N*sizeof(float));
st->in_mem = celt_alloc(N*C*sizeof(float));
@@ -130,7 +130,7 @@
ec_byte_writeclear(&st->buf);
mdct_clear(&st->mdct_lookup);
- spx_fft_destroy(st->fft);
+ free(st->fft);
celt_free(st->window);
celt_free(st->in_mem);
@@ -138,6 +138,10 @@
celt_free(st->out_mem);
celt_free(st->oldBandE);
+
+ celt_free(st->preemph_memE);
+ celt_free(st->preemph_memD);
+
alloc_clear(&st->alloc);
celt_free(st);
@@ -453,6 +457,9 @@
celt_free(st->out_mem);
celt_free(st->oldBandE);
+
+ celt_free(st->preemph_memD);
+
alloc_clear(&st->alloc);
celt_free(st);
--- a/libcelt/fftwrap.c
+++ /dev/null
@@ -1,299 +1,0 @@
-/* Copyright (C) 2005 Jean-Marc Valin
- File: fftwrap.c
-
- Wrapper for various FFTs
-
- Redistribution and use in source and binary forms, with or without
- modification, are permitted provided that the following conditions
- are met:
-
- - Redistributions of source code must retain the above copyright
- notice, this list of conditions and the following disclaimer.
-
- - Redistributions in binary form must reproduce the above copyright
- notice, this list of conditions and the following disclaimer in the
- documentation and/or other materials provided with the distribution.
-
- - Neither the name of the Xiph.org Foundation nor the names of its
- contributors may be used to endorse or promote products derived from
- this software without specific prior written permission.
-
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
- ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
- LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
- A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
- CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
-*/
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#define USE_SMALLFT
-/*#define USE_KISS_FFT*/
-
-
-#include "arch.h"
-#include "os_support.h"
-
-#ifdef FIXED_POINT
-static int maximize_range(spx_word16_t *in, spx_word16_t *out, spx_word16_t bound, int len)
-{- int i, shift;
- spx_word16_t max_val = 0;
- for (i=0;i<len;i++)
- {- if (in[i]>max_val)
- max_val = in[i];
- if (-in[i]>max_val)
- max_val = -in[i];
- }
- shift=0;
- while (max_val <= (bound>>1) && max_val != 0)
- {- max_val <<= 1;
- shift++;
- }
- for (i=0;i<len;i++)
- {- out[i] = in[i] << shift;
- }
- return shift;
-}
-
-static void renorm_range(spx_word16_t *in, spx_word16_t *out, int shift, int len)
-{- int i;
- for (i=0;i<len;i++)
- {- out[i] = (in[i] + (1<<(shift-1))) >> shift;
- }
-}
-#endif
-
-#ifdef USE_SMALLFT
-
-#include "smallft.h"
-#include <math.h>
-
-void *spx_fft_init(int size)
-{- struct drft_lookup *table;
- table = celt_alloc(sizeof(struct drft_lookup));
- spx_drft_init((struct drft_lookup *)table, size);
- return (void*)table;
-}
-
-void spx_fft_destroy(void *table)
-{- spx_drft_clear(table);
- celt_free(table);
-}
-
-void spx_fft(void *table, float *in, float *out)
-{- if (in==out)
- {- int i;
- celt_warning("FFT should not be done in-place");- float scale = 1./((struct drft_lookup *)table)->n;
- for (i=0;i<((struct drft_lookup *)table)->n;i++)
- out[i] = scale*in[i];
- } else {- int i;
- float scale = 1./((struct drft_lookup *)table)->n;
- for (i=0;i<((struct drft_lookup *)table)->n;i++)
- out[i] = scale*in[i];
- }
- spx_drft_forward((struct drft_lookup *)table, out);
-}
-
-void spx_ifft(void *table, float *in, float *out)
-{- if (in==out)
- {- celt_warning("FFT should not be done in-place");- } else {- int i;
- for (i=0;i<((struct drft_lookup *)table)->n;i++)
- out[i] = in[i];
- }
- spx_drft_backward((struct drft_lookup *)table, out);
-}
-
-#elif defined(USE_KISS_FFT)
-
-#include "kiss_fftr.h"
-#include "kiss_fft.h"
-
-struct kiss_config {- kiss_fftr_cfg forward;
- kiss_fftr_cfg backward;
- kiss_fft_cpx *freq_data;
- int N;
-};
-
-void *spx_fft_init(int size)
-{- struct kiss_config *table;
- table = speex_alloc(sizeof(struct kiss_config));
- table->freq_data = speex_alloc(sizeof(kiss_fft_cpx)*((size>>1)+1));
- table->forward = kiss_fftr_alloc(size,0,NULL,NULL);
- table->backward = kiss_fftr_alloc(size,1,NULL,NULL);
- table->N = size;
- return table;
-}
-
-void spx_fft_destroy(void *table)
-{- struct kiss_config *t = (struct kiss_config *)table;
- kiss_fftr_free(t->forward);
- kiss_fftr_free(t->backward);
- speex_free(t->freq_data);
- speex_free(table);
-}
-
-#ifdef FIXED_POINT
-
-void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
-{- int i;
- int shift;
- struct kiss_config *t = (struct kiss_config *)table;
- shift = maximize_range(in, in, 32000, t->N);
- kiss_fftr(t->forward, in, t->freq_data);
- out[0] = t->freq_data[0].r;
- for (i=1;i<t->N>>1;i++)
- {- out[(i<<1)-1] = t->freq_data[i].r;
- out[(i<<1)] = t->freq_data[i].i;
- }
- out[(i<<1)-1] = t->freq_data[i].r;
- renorm_range(in, in, shift, t->N);
- renorm_range(out, out, shift, t->N);
-}
-
-#else
-
-void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out)
-{- int i;
- float scale;
- struct kiss_config *t = (struct kiss_config *)table;
- scale = 1./t->N;
- kiss_fftr(t->forward, in, t->freq_data);
- out[0] = scale*t->freq_data[0].r;
- for (i=1;i<t->N>>1;i++)
- {- out[(i<<1)-1] = scale*t->freq_data[i].r;
- out[(i<<1)] = scale*t->freq_data[i].i;
- }
- out[(i<<1)-1] = scale*t->freq_data[i].r;
-}
-#endif
-
-void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out)
-{- int i;
- struct kiss_config *t = (struct kiss_config *)table;
- t->freq_data[0].r = in[0];
- t->freq_data[0].i = 0;
- for (i=1;i<t->N>>1;i++)
- {- t->freq_data[i].r = in[(i<<1)-1];
- t->freq_data[i].i = in[(i<<1)];
- }
- t->freq_data[i].r = in[(i<<1)-1];
- t->freq_data[i].i = 0;
-
- kiss_fftri(t->backward, t->freq_data, out);
-}
-
-
-#else
-
-#error No other FFT implemented
-
-#endif
-
-
-int fixed_point = 1;
-#ifdef FIXED_POINT
-#include "smallft.h"
-
-
-void spx_fft_float(void *table, float *in, float *out)
-{- int i;
-#ifdef USE_SMALLFT
- int N = ((struct drft_lookup *)table)->n;
-#elif defined(USE_KISS_FFT)
- int N = ((struct kiss_config *)table)->N;
-#else
-#endif
- spx_word16_t _in[N];
- spx_word16_t _out[N];
- for (i=0;i<N;i++)
- _in[i] = (int)floor(.5+in[i]);
- spx_fft(table, _in, _out);
- for (i=0;i<N;i++)
- out[i] = _out[i];
- if (!fixed_point)
- {- struct drft_lookup t;
- spx_drft_init(&t, ((struct kiss_config *)table)->N);
- float scale = 1./((struct kiss_config *)table)->N;
- for (i=0;i<((struct kiss_config *)table)->N;i++)
- out[i] = scale*in[i];
- spx_drft_forward(&t, out);
- spx_drft_clear(&t);
- }
-}
-
-void spx_ifft_float(void *table, float *in, float *out)
-{- int i;
-#ifdef USE_SMALLFT
- int N = ((struct drft_lookup *)table)->n;
-#elif defined(USE_KISS_FFT)
- int N = ((struct kiss_config *)table)->N;
-#else
-#endif
- spx_word16_t _in[N];
- spx_word16_t _out[N];
- for (i=0;i<N;i++)
- _in[i] = (int)floor(.5+in[i]);
- spx_ifft(table, _in, _out);
- for (i=0;i<N;i++)
- out[i] = _out[i];
- if (!fixed_point)
- {- int i;
- struct drft_lookup t;
- spx_drft_init(&t, ((struct kiss_config *)table)->N);
- for (i=0;i<((struct kiss_config *)table)->N;i++)
- out[i] = in[i];
- spx_drft_backward(&t, out);
- spx_drft_clear(&t);
- }
-}
-
-#else
-
-void spx_fft_float(void *table, float *in, float *out)
-{- spx_fft(table, in, out);
-}
-void spx_ifft_float(void *table, float *in, float *out)
-{- spx_ifft(table, in, out);
-}
-
-#endif
--- a/libcelt/fftwrap.h
+++ /dev/null
@@ -1,58 +1,0 @@
-/* Copyright (C) 2005 Jean-Marc Valin
- File: fftwrap.h
-
- Wrapper for various FFTs
-
- Redistribution and use in source and binary forms, with or without
- modification, are permitted provided that the following conditions
- are met:
-
- - Redistributions of source code must retain the above copyright
- notice, this list of conditions and the following disclaimer.
-
- - Redistributions in binary form must reproduce the above copyright
- notice, this list of conditions and the following disclaimer in the
- documentation and/or other materials provided with the distribution.
-
- - Neither the name of the Xiph.org Foundation nor the names of its
- contributors may be used to endorse or promote products derived from
- this software without specific prior written permission.
-
- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
- ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
- LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
- A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR
- CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
- EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
- PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
- PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
- LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
- NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
- SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
-
-*/
-
-#ifndef FFTWRAP_H
-#define FFTWRAP_H
-
-#include "arch.h"
-
-/** Compute tables for an FFT */
-void *spx_fft_init(int size);
-
-/** Destroy tables for an FFT */
-void spx_fft_destroy(void *table);
-
-/** Forward (real to half-complex) transform */
-void spx_fft(void *table, celt_word16_t *in, celt_word16_t *out);
-
-/** Backward (half-complex to real) transform */
-void spx_ifft(void *table, celt_word16_t *in, celt_word16_t *out);
-
-/** Forward (real to half-complex) transform of float data */
-void spx_fft_float(void *table, float *in, float *out);
-
-/** Backward (half-complex to real) transform of float data */
-void spx_ifft_float(void *table, float *in, float *out);
-
-#endif
--- a/libcelt/pitch.c
+++ b/libcelt/pitch.c
@@ -23,11 +23,10 @@
#include <stdio.h>
#include <math.h>
-#include "fftwrap.h"
#include "pitch.h"
#include "psy.h"
-void find_spectral_pitch(void *fft, float *x, float *y, int lag, int len, int C, int *pitch)
+void find_spectral_pitch(kiss_fftr_cfg fft, float *x, float *y, int lag, int len, int C, int *pitch)
{int c;
int n2 = lag/2;
@@ -48,12 +47,12 @@
yy[c*lag+i] = y[C*i+c];
}
- spx_fft(fft, xx, X);
- spx_fft(fft, yy, Y);
+ kiss_fftr(fft, xx, X);
+ kiss_fftr(fft, yy, Y);
+
compute_masking(X, curve, lag*C, 44100);
- X[0] = 0;
for (i=1;i<C*n2;i++)
{float n;
@@ -61,13 +60,12 @@
//n = 1;
n = 1.f/pow(1+curve[i],.5)/(i+60);
//n = 1.f/(1+curve[i]);
- float tmp = X[2*i-1];
- X[2*i-1] = (X[2*i-1]*Y[2*i-1] + X[2*i ]*Y[2*i ])*n;
- X[2*i ] = (- X[2*i ]*Y[2*i-1] + tmp*Y[2*i ])*n;
+ float tmp = X[2*i];
+ X[2*i] = (X[2*i ]*Y[2*i ] + X[2*i+1]*Y[2*i+1])*n;
+ X[2*i+1] = (- X[2*i+1]*Y[2*i ] + tmp*Y[2*i+1])*n;
}
- X[C*lag-1] = 0;
- X[0] = X[C*lag-1] = 0;
- spx_ifft(fft, X, xx);
+ X[0] = X[1] = 0;
+ kiss_fftri(fft, X, xx);
float max_corr=-1e10;
//int pitch;
--- a/libcelt/pitch.h
+++ b/libcelt/pitch.h
@@ -23,6 +23,8 @@
#ifndef _PITCH_H
#define _PITCH_H
-void find_spectral_pitch(void *fft, float *x, float *y, int lag, int len, int C, int *pitch);
+#include "kiss_fftr.h"
+
+void find_spectral_pitch(kiss_fftr_cfg fft, float *x, float *y, int lag, int len, int C, int *pitch);
#endif
--- a/libcelt/psy.c
+++ b/libcelt/psy.c
@@ -111,7 +111,7 @@
float psd[N];
psd[0] = X[0]*X[0];
for (i=1;i<N;i++)
- psd[i] = X[i*2-1]*X[i*2-1] + X[i*2]*X[i*2];
+ psd[i] = X[i*2]*X[i*2] + X[i*2+1]*X[i*2+1];
/* TODO: Do tone masking */
/* Noise masking */
spreading_func(psd, mask, N, Fs);
--- a/libcelt/smallft.c
+++ /dev/null
@@ -1,1260 +1,0 @@
-/********************************************************************
- * *
- * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
- * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
- * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
- * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
- * *
- * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 *
- * by the XIPHOPHORUS Company http://www.xiph.org/ *
- * *
- ********************************************************************
-
- function: *unnormalized* fft transform
- last mod: $Id: smallft.c,v 1.19 2003/10/08 05:12:37 jm Exp $
-
- ********************************************************************/
-
-/* FFT implementation from OggSquish, minus cosine transforms,
- * minus all but radix 2/4 case. In Vorbis we only need this
- * cut-down version.
- *
- * To do more than just power-of-two sized vectors, see the full
- * version I wrote for NetLib.
- *
- * Note that the packing is a little strange; rather than the FFT r/i
- * packing following R_0, I_n, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1,
- * it follows R_0, R_1, I_1, R_2, I_2 ... R_n-1, I_n-1, I_n like the
- * FORTRAN version
- */
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <math.h>
-#include "smallft.h"
-#include "os_support.h"
-
-static void drfti1(int n, float *wa, int *ifac){- static int ntryh[4] = { 4,2,3,5 };- static float tpi = 6.28318530717958648f;
- float arg,argh,argld,fi;
- int ntry=0,i,j=-1;
- int k1, l1, l2, ib;
- int ld, ii, ip, is, nq, nr;
- int ido, ipm, nfm1;
- int nl=n;
- int nf=0;
-
- L101:
- j++;
- if (j < 4)
- ntry=ntryh[j];
- else
- ntry+=2;
-
- L104:
- nq=nl/ntry;
- nr=nl-ntry*nq;
- if (nr!=0) goto L101;
-
- nf++;
- ifac[nf+1]=ntry;
- nl=nq;
- if(ntry!=2)goto L107;
- if(nf==1)goto L107;
-
- for (i=1;i<nf;i++){- ib=nf-i+1;
- ifac[ib+1]=ifac[ib];
- }
- ifac[2] = 2;
-
- L107:
- if(nl!=1)goto L104;
- ifac[0]=n;
- ifac[1]=nf;
- argh=tpi/n;
- is=0;
- nfm1=nf-1;
- l1=1;
-
- if(nfm1==0)return;
-
- for (k1=0;k1<nfm1;k1++){- ip=ifac[k1+2];
- ld=0;
- l2=l1*ip;
- ido=n/l2;
- ipm=ip-1;
-
- for (j=0;j<ipm;j++){- ld+=l1;
- i=is;
- argld=(float)ld*argh;
- fi=0.f;
- for (ii=2;ii<ido;ii+=2){- fi+=1.f;
- arg=fi*argld;
- wa[i++]=cos(arg);
- wa[i++]=sin(arg);
- }
- is+=ido;
- }
- l1=l2;
- }
-}
-
-static void fdrffti(int n, float *wsave, int *ifac){-
- if (n == 1) return;
- drfti1(n, wsave+n, ifac);
-}
-
-static void dradf2(int ido,int l1,float *cc,float *ch,float *wa1){- int i,k;
- float ti2,tr2;
- int t0,t1,t2,t3,t4,t5,t6;
-
- t1=0;
- t0=(t2=l1*ido);
- t3=ido<<1;
- for(k=0;k<l1;k++){- ch[t1<<1]=cc[t1]+cc[t2];
- ch[(t1<<1)+t3-1]=cc[t1]-cc[t2];
- t1+=ido;
- t2+=ido;
- }
-
- if(ido<2)return;
- if(ido==2)goto L105;
-
- t1=0;
- t2=t0;
- for(k=0;k<l1;k++){- t3=t2;
- t4=(t1<<1)+(ido<<1);
- t5=t1;
- t6=t1+t1;
- for(i=2;i<ido;i+=2){- t3+=2;
- t4-=2;
- t5+=2;
- t6+=2;
- tr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
- ti2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
- ch[t6]=cc[t5]+ti2;
- ch[t4]=ti2-cc[t5];
- ch[t6-1]=cc[t5-1]+tr2;
- ch[t4-1]=cc[t5-1]-tr2;
- }
- t1+=ido;
- t2+=ido;
- }
-
- if(ido%2==1)return;
-
- L105:
- t3=(t2=(t1=ido)-1);
- t2+=t0;
- for(k=0;k<l1;k++){- ch[t1]=-cc[t2];
- ch[t1-1]=cc[t3];
- t1+=ido<<1;
- t2+=ido;
- t3+=ido;
- }
-}
-
-static void dradf4(int ido,int l1,float *cc,float *ch,float *wa1,
- float *wa2,float *wa3){- static float hsqt2 = .70710678118654752f;
- int i,k,t0,t1,t2,t3,t4,t5,t6;
- float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
- t0=l1*ido;
-
- t1=t0;
- t4=t1<<1;
- t2=t1+(t1<<1);
- t3=0;
-
- for(k=0;k<l1;k++){- tr1=cc[t1]+cc[t2];
- tr2=cc[t3]+cc[t4];
-
- ch[t5=t3<<2]=tr1+tr2;
- ch[(ido<<2)+t5-1]=tr2-tr1;
- ch[(t5+=(ido<<1))-1]=cc[t3]-cc[t4];
- ch[t5]=cc[t2]-cc[t1];
-
- t1+=ido;
- t2+=ido;
- t3+=ido;
- t4+=ido;
- }
-
- if(ido<2)return;
- if(ido==2)goto L105;
-
-
- t1=0;
- for(k=0;k<l1;k++){- t2=t1;
- t4=t1<<2;
- t5=(t6=ido<<1)+t4;
- for(i=2;i<ido;i+=2){- t3=(t2+=2);
- t4+=2;
- t5-=2;
-
- t3+=t0;
- cr2=wa1[i-2]*cc[t3-1]+wa1[i-1]*cc[t3];
- ci2=wa1[i-2]*cc[t3]-wa1[i-1]*cc[t3-1];
- t3+=t0;
- cr3=wa2[i-2]*cc[t3-1]+wa2[i-1]*cc[t3];
- ci3=wa2[i-2]*cc[t3]-wa2[i-1]*cc[t3-1];
- t3+=t0;
- cr4=wa3[i-2]*cc[t3-1]+wa3[i-1]*cc[t3];
- ci4=wa3[i-2]*cc[t3]-wa3[i-1]*cc[t3-1];
-
- tr1=cr2+cr4;
- tr4=cr4-cr2;
- ti1=ci2+ci4;
- ti4=ci2-ci4;
-
- ti2=cc[t2]+ci3;
- ti3=cc[t2]-ci3;
- tr2=cc[t2-1]+cr3;
- tr3=cc[t2-1]-cr3;
-
- ch[t4-1]=tr1+tr2;
- ch[t4]=ti1+ti2;
-
- ch[t5-1]=tr3-ti4;
- ch[t5]=tr4-ti3;
-
- ch[t4+t6-1]=ti4+tr3;
- ch[t4+t6]=tr4+ti3;
-
- ch[t5+t6-1]=tr2-tr1;
- ch[t5+t6]=ti1-ti2;
- }
- t1+=ido;
- }
- if(ido&1)return;
-
- L105:
-
- t2=(t1=t0+ido-1)+(t0<<1);
- t3=ido<<2;
- t4=ido;
- t5=ido<<1;
- t6=ido;
-
- for(k=0;k<l1;k++){- ti1=-hsqt2*(cc[t1]+cc[t2]);
- tr1=hsqt2*(cc[t1]-cc[t2]);
-
- ch[t4-1]=tr1+cc[t6-1];
- ch[t4+t5-1]=cc[t6-1]-tr1;
-
- ch[t4]=ti1-cc[t1+t0];
- ch[t4+t5]=ti1+cc[t1+t0];
-
- t1+=ido;
- t2+=ido;
- t4+=t3;
- t6+=ido;
- }
-}
-
-static void dradfg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
- float *c2,float *ch,float *ch2,float *wa){-
- static float tpi=6.283185307179586f;
- int idij,ipph,i,j,k,l,ic,ik,is;
- int t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
- float dc2,ai1,ai2,ar1,ar2,ds2;
- int nbd;
- float dcp,arg,dsp,ar1h,ar2h;
- int idp2,ipp2;
-
- arg=tpi/(float)ip;
- dcp=cos(arg);
- dsp=sin(arg);
- ipph=(ip+1)>>1;
- ipp2=ip;
- idp2=ido;
- nbd=(ido-1)>>1;
- t0=l1*ido;
- t10=ip*ido;
-
- if(ido==1)goto L119;
- for(ik=0;ik<idl1;ik++)ch2[ik]=c2[ik];
-
- t1=0;
- for(j=1;j<ip;j++){- t1+=t0;
- t2=t1;
- for(k=0;k<l1;k++){- ch[t2]=c1[t2];
- t2+=ido;
- }
- }
-
- is=-ido;
- t1=0;
- if(nbd>l1){- for(j=1;j<ip;j++){- t1+=t0;
- is+=ido;
- t2= -ido+t1;
- for(k=0;k<l1;k++){- idij=is-1;
- t2+=ido;
- t3=t2;
- for(i=2;i<ido;i+=2){- idij+=2;
- t3+=2;
- ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
- ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
- }
- }
- }
- }else{-
- for(j=1;j<ip;j++){- is+=ido;
- idij=is-1;
- t1+=t0;
- t2=t1;
- for(i=2;i<ido;i+=2){- idij+=2;
- t2+=2;
- t3=t2;
- for(k=0;k<l1;k++){- ch[t3-1]=wa[idij-1]*c1[t3-1]+wa[idij]*c1[t3];
- ch[t3]=wa[idij-1]*c1[t3]-wa[idij]*c1[t3-1];
- t3+=ido;
- }
- }
- }
- }
-
- t1=0;
- t2=ipp2*t0;
- if(nbd<l1){- for(j=1;j<ipph;j++){- t1+=t0;
- t2-=t0;
- t3=t1;
- t4=t2;
- for(i=2;i<ido;i+=2){- t3+=2;
- t4+=2;
- t5=t3-ido;
- t6=t4-ido;
- for(k=0;k<l1;k++){- t5+=ido;
- t6+=ido;
- c1[t5-1]=ch[t5-1]+ch[t6-1];
- c1[t6-1]=ch[t5]-ch[t6];
- c1[t5]=ch[t5]+ch[t6];
- c1[t6]=ch[t6-1]-ch[t5-1];
- }
- }
- }
- }else{- for(j=1;j<ipph;j++){- t1+=t0;
- t2-=t0;
- t3=t1;
- t4=t2;
- for(k=0;k<l1;k++){- t5=t3;
- t6=t4;
- for(i=2;i<ido;i+=2){- t5+=2;
- t6+=2;
- c1[t5-1]=ch[t5-1]+ch[t6-1];
- c1[t6-1]=ch[t5]-ch[t6];
- c1[t5]=ch[t5]+ch[t6];
- c1[t6]=ch[t6-1]-ch[t5-1];
- }
- t3+=ido;
- t4+=ido;
- }
- }
- }
-
-L119:
- for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
-
- t1=0;
- t2=ipp2*idl1;
- for(j=1;j<ipph;j++){- t1+=t0;
- t2-=t0;
- t3=t1-ido;
- t4=t2-ido;
- for(k=0;k<l1;k++){- t3+=ido;
- t4+=ido;
- c1[t3]=ch[t3]+ch[t4];
- c1[t4]=ch[t4]-ch[t3];
- }
- }
-
- ar1=1.f;
- ai1=0.f;
- t1=0;
- t2=ipp2*idl1;
- t3=(ip-1)*idl1;
- for(l=1;l<ipph;l++){- t1+=idl1;
- t2-=idl1;
- ar1h=dcp*ar1-dsp*ai1;
- ai1=dcp*ai1+dsp*ar1;
- ar1=ar1h;
- t4=t1;
- t5=t2;
- t6=t3;
- t7=idl1;
-
- for(ik=0;ik<idl1;ik++){- ch2[t4++]=c2[ik]+ar1*c2[t7++];
- ch2[t5++]=ai1*c2[t6++];
- }
-
- dc2=ar1;
- ds2=ai1;
- ar2=ar1;
- ai2=ai1;
-
- t4=idl1;
- t5=(ipp2-1)*idl1;
- for(j=2;j<ipph;j++){- t4+=idl1;
- t5-=idl1;
-
- ar2h=dc2*ar2-ds2*ai2;
- ai2=dc2*ai2+ds2*ar2;
- ar2=ar2h;
-
- t6=t1;
- t7=t2;
- t8=t4;
- t9=t5;
- for(ik=0;ik<idl1;ik++){- ch2[t6++]+=ar2*c2[t8++];
- ch2[t7++]+=ai2*c2[t9++];
- }
- }
- }
-
- t1=0;
- for(j=1;j<ipph;j++){- t1+=idl1;
- t2=t1;
- for(ik=0;ik<idl1;ik++)ch2[ik]+=c2[t2++];
- }
-
- if(ido<l1)goto L132;
-
- t1=0;
- t2=0;
- for(k=0;k<l1;k++){- t3=t1;
- t4=t2;
- for(i=0;i<ido;i++)cc[t4++]=ch[t3++];
- t1+=ido;
- t2+=t10;
- }
-
- goto L135;
-
- L132:
- for(i=0;i<ido;i++){- t1=i;
- t2=i;
- for(k=0;k<l1;k++){- cc[t2]=ch[t1];
- t1+=ido;
- t2+=t10;
- }
- }
-
- L135:
- t1=0;
- t2=ido<<1;
- t3=0;
- t4=ipp2*t0;
- for(j=1;j<ipph;j++){-
- t1+=t2;
- t3+=t0;
- t4-=t0;
-
- t5=t1;
- t6=t3;
- t7=t4;
-
- for(k=0;k<l1;k++){- cc[t5-1]=ch[t6];
- cc[t5]=ch[t7];
- t5+=t10;
- t6+=ido;
- t7+=ido;
- }
- }
-
- if(ido==1)return;
- if(nbd<l1)goto L141;
-
- t1=-ido;
- t3=0;
- t4=0;
- t5=ipp2*t0;
- for(j=1;j<ipph;j++){- t1+=t2;
- t3+=t2;
- t4+=t0;
- t5-=t0;
- t6=t1;
- t7=t3;
- t8=t4;
- t9=t5;
- for(k=0;k<l1;k++){- for(i=2;i<ido;i+=2){- ic=idp2-i;
- cc[i+t7-1]=ch[i+t8-1]+ch[i+t9-1];
- cc[ic+t6-1]=ch[i+t8-1]-ch[i+t9-1];
- cc[i+t7]=ch[i+t8]+ch[i+t9];
- cc[ic+t6]=ch[i+t9]-ch[i+t8];
- }
- t6+=t10;
- t7+=t10;
- t8+=ido;
- t9+=ido;
- }
- }
- return;
-
- L141:
-
- t1=-ido;
- t3=0;
- t4=0;
- t5=ipp2*t0;
- for(j=1;j<ipph;j++){- t1+=t2;
- t3+=t2;
- t4+=t0;
- t5-=t0;
- for(i=2;i<ido;i+=2){- t6=idp2+t1-i;
- t7=i+t3;
- t8=i+t4;
- t9=i+t5;
- for(k=0;k<l1;k++){- cc[t7-1]=ch[t8-1]+ch[t9-1];
- cc[t6-1]=ch[t8-1]-ch[t9-1];
- cc[t7]=ch[t8]+ch[t9];
- cc[t6]=ch[t9]-ch[t8];
- t6+=t10;
- t7+=t10;
- t8+=ido;
- t9+=ido;
- }
- }
- }
-}
-
-static void drftf1(int n,float *c,float *ch,float *wa,int *ifac){- int i,k1,l1,l2;
- int na,kh,nf;
- int ip,iw,ido,idl1,ix2,ix3;
-
- nf=ifac[1];
- na=1;
- l2=n;
- iw=n;
-
- for(k1=0;k1<nf;k1++){- kh=nf-k1;
- ip=ifac[kh+1];
- l1=l2/ip;
- ido=n/l2;
- idl1=ido*l1;
- iw-=(ip-1)*ido;
- na=1-na;
-
- if(ip!=4)goto L102;
-
- ix2=iw+ido;
- ix3=ix2+ido;
- if(na!=0)
- dradf4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
- else
- dradf4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
- goto L110;
-
- L102:
- if(ip!=2)goto L104;
- if(na!=0)goto L103;
-
- dradf2(ido,l1,c,ch,wa+iw-1);
- goto L110;
-
- L103:
- dradf2(ido,l1,ch,c,wa+iw-1);
- goto L110;
-
- L104:
- if(ido==1)na=1-na;
- if(na!=0)goto L109;
-
- dradfg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
- na=1;
- goto L110;
-
- L109:
- dradfg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
- na=0;
-
- L110:
- l2=l1;
- }
-
- if(na==1)return;
-
- for(i=0;i<n;i++)c[i]=ch[i];
-}
-
-static void dradb2(int ido,int l1,float *cc,float *ch,float *wa1){- int i,k,t0,t1,t2,t3,t4,t5,t6;
- float ti2,tr2;
-
- t0=l1*ido;
-
- t1=0;
- t2=0;
- t3=(ido<<1)-1;
- for(k=0;k<l1;k++){- ch[t1]=cc[t2]+cc[t3+t2];
- ch[t1+t0]=cc[t2]-cc[t3+t2];
- t2=(t1+=ido)<<1;
- }
-
- if(ido<2)return;
- if(ido==2)goto L105;
-
- t1=0;
- t2=0;
- for(k=0;k<l1;k++){- t3=t1;
- t5=(t4=t2)+(ido<<1);
- t6=t0+t1;
- for(i=2;i<ido;i+=2){- t3+=2;
- t4+=2;
- t5-=2;
- t6+=2;
- ch[t3-1]=cc[t4-1]+cc[t5-1];
- tr2=cc[t4-1]-cc[t5-1];
- ch[t3]=cc[t4]-cc[t5];
- ti2=cc[t4]+cc[t5];
- ch[t6-1]=wa1[i-2]*tr2-wa1[i-1]*ti2;
- ch[t6]=wa1[i-2]*ti2+wa1[i-1]*tr2;
- }
- t2=(t1+=ido)<<1;
- }
-
- if(ido%2==1)return;
-
-L105:
- t1=ido-1;
- t2=ido-1;
- for(k=0;k<l1;k++){- ch[t1]=cc[t2]+cc[t2];
- ch[t1+t0]=-(cc[t2+1]+cc[t2+1]);
- t1+=ido;
- t2+=ido<<1;
- }
-}
-
-static void dradb3(int ido,int l1,float *cc,float *ch,float *wa1,
- float *wa2){- static float taur = -.5f;
- static float taui = .8660254037844386f;
- int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10;
- float ci2,ci3,di2,di3,cr2,cr3,dr2,dr3,ti2,tr2;
- t0=l1*ido;
-
- t1=0;
- t2=t0<<1;
- t3=ido<<1;
- t4=ido+(ido<<1);
- t5=0;
- for(k=0;k<l1;k++){- tr2=cc[t3-1]+cc[t3-1];
- cr2=cc[t5]+(taur*tr2);
- ch[t1]=cc[t5]+tr2;
- ci3=taui*(cc[t3]+cc[t3]);
- ch[t1+t0]=cr2-ci3;
- ch[t1+t2]=cr2+ci3;
- t1+=ido;
- t3+=t4;
- t5+=t4;
- }
-
- if(ido==1)return;
-
- t1=0;
- t3=ido<<1;
- for(k=0;k<l1;k++){- t7=t1+(t1<<1);
- t6=(t5=t7+t3);
- t8=t1;
- t10=(t9=t1+t0)+t0;
-
- for(i=2;i<ido;i+=2){- t5+=2;
- t6-=2;
- t7+=2;
- t8+=2;
- t9+=2;
- t10+=2;
- tr2=cc[t5-1]+cc[t6-1];
- cr2=cc[t7-1]+(taur*tr2);
- ch[t8-1]=cc[t7-1]+tr2;
- ti2=cc[t5]-cc[t6];
- ci2=cc[t7]+(taur*ti2);
- ch[t8]=cc[t7]+ti2;
- cr3=taui*(cc[t5-1]-cc[t6-1]);
- ci3=taui*(cc[t5]+cc[t6]);
- dr2=cr2-ci3;
- dr3=cr2+ci3;
- di2=ci2+cr3;
- di3=ci2-cr3;
- ch[t9-1]=wa1[i-2]*dr2-wa1[i-1]*di2;
- ch[t9]=wa1[i-2]*di2+wa1[i-1]*dr2;
- ch[t10-1]=wa2[i-2]*dr3-wa2[i-1]*di3;
- ch[t10]=wa2[i-2]*di3+wa2[i-1]*dr3;
- }
- t1+=ido;
- }
-}
-
-static void dradb4(int ido,int l1,float *cc,float *ch,float *wa1,
- float *wa2,float *wa3){- static float sqrt2=1.414213562373095f;
- int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
- float ci2,ci3,ci4,cr2,cr3,cr4,ti1,ti2,ti3,ti4,tr1,tr2,tr3,tr4;
- t0=l1*ido;
-
- t1=0;
- t2=ido<<2;
- t3=0;
- t6=ido<<1;
- for(k=0;k<l1;k++){- t4=t3+t6;
- t5=t1;
- tr3=cc[t4-1]+cc[t4-1];
- tr4=cc[t4]+cc[t4];
- tr1=cc[t3]-cc[(t4+=t6)-1];
- tr2=cc[t3]+cc[t4-1];
- ch[t5]=tr2+tr3;
- ch[t5+=t0]=tr1-tr4;
- ch[t5+=t0]=tr2-tr3;
- ch[t5+=t0]=tr1+tr4;
- t1+=ido;
- t3+=t2;
- }
-
- if(ido<2)return;
- if(ido==2)goto L105;
-
- t1=0;
- for(k=0;k<l1;k++){- t5=(t4=(t3=(t2=t1<<2)+t6))+t6;
- t7=t1;
- for(i=2;i<ido;i+=2){- t2+=2;
- t3+=2;
- t4-=2;
- t5-=2;
- t7+=2;
- ti1=cc[t2]+cc[t5];
- ti2=cc[t2]-cc[t5];
- ti3=cc[t3]-cc[t4];
- tr4=cc[t3]+cc[t4];
- tr1=cc[t2-1]-cc[t5-1];
- tr2=cc[t2-1]+cc[t5-1];
- ti4=cc[t3-1]-cc[t4-1];
- tr3=cc[t3-1]+cc[t4-1];
- ch[t7-1]=tr2+tr3;
- cr3=tr2-tr3;
- ch[t7]=ti2+ti3;
- ci3=ti2-ti3;
- cr2=tr1-tr4;
- cr4=tr1+tr4;
- ci2=ti1+ti4;
- ci4=ti1-ti4;
-
- ch[(t8=t7+t0)-1]=wa1[i-2]*cr2-wa1[i-1]*ci2;
- ch[t8]=wa1[i-2]*ci2+wa1[i-1]*cr2;
- ch[(t8+=t0)-1]=wa2[i-2]*cr3-wa2[i-1]*ci3;
- ch[t8]=wa2[i-2]*ci3+wa2[i-1]*cr3;
- ch[(t8+=t0)-1]=wa3[i-2]*cr4-wa3[i-1]*ci4;
- ch[t8]=wa3[i-2]*ci4+wa3[i-1]*cr4;
- }
- t1+=ido;
- }
-
- if(ido%2 == 1)return;
-
- L105:
-
- t1=ido;
- t2=ido<<2;
- t3=ido-1;
- t4=ido+(ido<<1);
- for(k=0;k<l1;k++){- t5=t3;
- ti1=cc[t1]+cc[t4];
- ti2=cc[t4]-cc[t1];
- tr1=cc[t1-1]-cc[t4-1];
- tr2=cc[t1-1]+cc[t4-1];
- ch[t5]=tr2+tr2;
- ch[t5+=t0]=sqrt2*(tr1-ti1);
- ch[t5+=t0]=ti2+ti2;
- ch[t5+=t0]=-sqrt2*(tr1+ti1);
-
- t3+=ido;
- t1+=t2;
- t4+=t2;
- }
-}
-
-static void dradbg(int ido,int ip,int l1,int idl1,float *cc,float *c1,
- float *c2,float *ch,float *ch2,float *wa){- static float tpi=6.283185307179586f;
- int idij,ipph,i,j,k,l,ik,is,t0,t1,t2,t3,t4,t5,t6,t7,t8,t9,t10,
- t11,t12;
- float dc2,ai1,ai2,ar1,ar2,ds2;
- int nbd;
- float dcp,arg,dsp,ar1h,ar2h;
- int ipp2;
-
- t10=ip*ido;
- t0=l1*ido;
- arg=tpi/(float)ip;
- dcp=cos(arg);
- dsp=sin(arg);
- nbd=(ido-1)>>1;
- ipp2=ip;
- ipph=(ip+1)>>1;
- if(ido<l1)goto L103;
-
- t1=0;
- t2=0;
- for(k=0;k<l1;k++){- t3=t1;
- t4=t2;
- for(i=0;i<ido;i++){- ch[t3]=cc[t4];
- t3++;
- t4++;
- }
- t1+=ido;
- t2+=t10;
- }
- goto L106;
-
- L103:
- t1=0;
- for(i=0;i<ido;i++){- t2=t1;
- t3=t1;
- for(k=0;k<l1;k++){- ch[t2]=cc[t3];
- t2+=ido;
- t3+=t10;
- }
- t1++;
- }
-
- L106:
- t1=0;
- t2=ipp2*t0;
- t7=(t5=ido<<1);
- for(j=1;j<ipph;j++){- t1+=t0;
- t2-=t0;
- t3=t1;
- t4=t2;
- t6=t5;
- for(k=0;k<l1;k++){- ch[t3]=cc[t6-1]+cc[t6-1];
- ch[t4]=cc[t6]+cc[t6];
- t3+=ido;
- t4+=ido;
- t6+=t10;
- }
- t5+=t7;
- }
-
- if (ido == 1)goto L116;
- if(nbd<l1)goto L112;
-
- t1=0;
- t2=ipp2*t0;
- t7=0;
- for(j=1;j<ipph;j++){- t1+=t0;
- t2-=t0;
- t3=t1;
- t4=t2;
-
- t7+=(ido<<1);
- t8=t7;
- for(k=0;k<l1;k++){- t5=t3;
- t6=t4;
- t9=t8;
- t11=t8;
- for(i=2;i<ido;i+=2){- t5+=2;
- t6+=2;
- t9+=2;
- t11-=2;
- ch[t5-1]=cc[t9-1]+cc[t11-1];
- ch[t6-1]=cc[t9-1]-cc[t11-1];
- ch[t5]=cc[t9]-cc[t11];
- ch[t6]=cc[t9]+cc[t11];
- }
- t3+=ido;
- t4+=ido;
- t8+=t10;
- }
- }
- goto L116;
-
- L112:
- t1=0;
- t2=ipp2*t0;
- t7=0;
- for(j=1;j<ipph;j++){- t1+=t0;
- t2-=t0;
- t3=t1;
- t4=t2;
- t7+=(ido<<1);
- t8=t7;
- t9=t7;
- for(i=2;i<ido;i+=2){- t3+=2;
- t4+=2;
- t8+=2;
- t9-=2;
- t5=t3;
- t6=t4;
- t11=t8;
- t12=t9;
- for(k=0;k<l1;k++){- ch[t5-1]=cc[t11-1]+cc[t12-1];
- ch[t6-1]=cc[t11-1]-cc[t12-1];
- ch[t5]=cc[t11]-cc[t12];
- ch[t6]=cc[t11]+cc[t12];
- t5+=ido;
- t6+=ido;
- t11+=t10;
- t12+=t10;
- }
- }
- }
-
-L116:
- ar1=1.f;
- ai1=0.f;
- t1=0;
- t9=(t2=ipp2*idl1);
- t3=(ip-1)*idl1;
- for(l=1;l<ipph;l++){- t1+=idl1;
- t2-=idl1;
-
- ar1h=dcp*ar1-dsp*ai1;
- ai1=dcp*ai1+dsp*ar1;
- ar1=ar1h;
- t4=t1;
- t5=t2;
- t6=0;
- t7=idl1;
- t8=t3;
- for(ik=0;ik<idl1;ik++){- c2[t4++]=ch2[t6++]+ar1*ch2[t7++];
- c2[t5++]=ai1*ch2[t8++];
- }
- dc2=ar1;
- ds2=ai1;
- ar2=ar1;
- ai2=ai1;
-
- t6=idl1;
- t7=t9-idl1;
- for(j=2;j<ipph;j++){- t6+=idl1;
- t7-=idl1;
- ar2h=dc2*ar2-ds2*ai2;
- ai2=dc2*ai2+ds2*ar2;
- ar2=ar2h;
- t4=t1;
- t5=t2;
- t11=t6;
- t12=t7;
- for(ik=0;ik<idl1;ik++){- c2[t4++]+=ar2*ch2[t11++];
- c2[t5++]+=ai2*ch2[t12++];
- }
- }
- }
-
- t1=0;
- for(j=1;j<ipph;j++){- t1+=idl1;
- t2=t1;
- for(ik=0;ik<idl1;ik++)ch2[ik]+=ch2[t2++];
- }
-
- t1=0;
- t2=ipp2*t0;
- for(j=1;j<ipph;j++){- t1+=t0;
- t2-=t0;
- t3=t1;
- t4=t2;
- for(k=0;k<l1;k++){- ch[t3]=c1[t3]-c1[t4];
- ch[t4]=c1[t3]+c1[t4];
- t3+=ido;
- t4+=ido;
- }
- }
-
- if(ido==1)goto L132;
- if(nbd<l1)goto L128;
-
- t1=0;
- t2=ipp2*t0;
- for(j=1;j<ipph;j++){- t1+=t0;
- t2-=t0;
- t3=t1;
- t4=t2;
- for(k=0;k<l1;k++){- t5=t3;
- t6=t4;
- for(i=2;i<ido;i+=2){- t5+=2;
- t6+=2;
- ch[t5-1]=c1[t5-1]-c1[t6];
- ch[t6-1]=c1[t5-1]+c1[t6];
- ch[t5]=c1[t5]+c1[t6-1];
- ch[t6]=c1[t5]-c1[t6-1];
- }
- t3+=ido;
- t4+=ido;
- }
- }
- goto L132;
-
- L128:
- t1=0;
- t2=ipp2*t0;
- for(j=1;j<ipph;j++){- t1+=t0;
- t2-=t0;
- t3=t1;
- t4=t2;
- for(i=2;i<ido;i+=2){- t3+=2;
- t4+=2;
- t5=t3;
- t6=t4;
- for(k=0;k<l1;k++){- ch[t5-1]=c1[t5-1]-c1[t6];
- ch[t6-1]=c1[t5-1]+c1[t6];
- ch[t5]=c1[t5]+c1[t6-1];
- ch[t6]=c1[t5]-c1[t6-1];
- t5+=ido;
- t6+=ido;
- }
- }
- }
-
-L132:
- if(ido==1)return;
-
- for(ik=0;ik<idl1;ik++)c2[ik]=ch2[ik];
-
- t1=0;
- for(j=1;j<ip;j++){- t2=(t1+=t0);
- for(k=0;k<l1;k++){- c1[t2]=ch[t2];
- t2+=ido;
- }
- }
-
- if(nbd>l1)goto L139;
-
- is= -ido-1;
- t1=0;
- for(j=1;j<ip;j++){- is+=ido;
- t1+=t0;
- idij=is;
- t2=t1;
- for(i=2;i<ido;i+=2){- t2+=2;
- idij+=2;
- t3=t2;
- for(k=0;k<l1;k++){- c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
- c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
- t3+=ido;
- }
- }
- }
- return;
-
- L139:
- is= -ido-1;
- t1=0;
- for(j=1;j<ip;j++){- is+=ido;
- t1+=t0;
- t2=t1;
- for(k=0;k<l1;k++){- idij=is;
- t3=t2;
- for(i=2;i<ido;i+=2){- idij+=2;
- t3+=2;
- c1[t3-1]=wa[idij-1]*ch[t3-1]-wa[idij]*ch[t3];
- c1[t3]=wa[idij-1]*ch[t3]+wa[idij]*ch[t3-1];
- }
- t2+=ido;
- }
- }
-}
-
-static void drftb1(int n, float *c, float *ch, float *wa, int *ifac){- int i,k1,l1,l2;
- int na;
- int nf,ip,iw,ix2,ix3,ido,idl1;
-
- nf=ifac[1];
- na=0;
- l1=1;
- iw=1;
-
- for(k1=0;k1<nf;k1++){- ip=ifac[k1 + 2];
- l2=ip*l1;
- ido=n/l2;
- idl1=ido*l1;
- if(ip!=4)goto L103;
- ix2=iw+ido;
- ix3=ix2+ido;
-
- if(na!=0)
- dradb4(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1);
- else
- dradb4(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1);
- na=1-na;
- goto L115;
-
- L103:
- if(ip!=2)goto L106;
-
- if(na!=0)
- dradb2(ido,l1,ch,c,wa+iw-1);
- else
- dradb2(ido,l1,c,ch,wa+iw-1);
- na=1-na;
- goto L115;
-
- L106:
- if(ip!=3)goto L109;
-
- ix2=iw+ido;
- if(na!=0)
- dradb3(ido,l1,ch,c,wa+iw-1,wa+ix2-1);
- else
- dradb3(ido,l1,c,ch,wa+iw-1,wa+ix2-1);
- na=1-na;
- goto L115;
-
- L109:
-/* The radix five case can be translated later..... */
-/* if(ip!=5)goto L112;
-
- ix2=iw+ido;
- ix3=ix2+ido;
- ix4=ix3+ido;
- if(na!=0)
- dradb5(ido,l1,ch,c,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
- else
- dradb5(ido,l1,c,ch,wa+iw-1,wa+ix2-1,wa+ix3-1,wa+ix4-1);
- na=1-na;
- goto L115;
-
- L112:*/
- if(na!=0)
- dradbg(ido,ip,l1,idl1,ch,ch,ch,c,c,wa+iw-1);
- else
- dradbg(ido,ip,l1,idl1,c,c,c,ch,ch,wa+iw-1);
- if(ido==1)na=1-na;
-
- L115:
- l1=l2;
- iw+=(ip-1)*ido;
- }
-
- if(na==0)return;
-
- for(i=0;i<n;i++)c[i]=ch[i];
-}
-
-void spx_drft_forward(struct drft_lookup *l,float *data){- if(l->n==1)return;
- drftf1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
-}
-
-void spx_drft_backward(struct drft_lookup *l,float *data){- if (l->n==1)return;
- drftb1(l->n,data,l->trigcache,l->trigcache+l->n,l->splitcache);
-}
-
-void spx_drft_init(struct drft_lookup *l,int n)
-{- l->n=n;
- l->trigcache=(float*)celt_alloc(3*n*sizeof(*l->trigcache));
- l->splitcache=(int*)celt_alloc(32*sizeof(*l->splitcache));
- fdrffti(n, l->trigcache, l->splitcache);
-}
-
-void spx_drft_clear(struct drft_lookup *l)
-{- if(l)
- {- if(l->trigcache)
- celt_free(l->trigcache);
- if(l->splitcache)
- celt_free(l->splitcache);
- }
-}
--- a/libcelt/smallft.h
+++ /dev/null
@@ -1,46 +1,0 @@
-/********************************************************************
- * *
- * THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. *
- * USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS *
- * GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE *
- * IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. *
- * *
- * THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2001 *
- * by the XIPHOPHORUS Company http://www.xiph.org/ *
- * *
- ********************************************************************
-
- function: fft transform
- last mod: $Id: smallft.h,v 1.3 2003/09/16 18:35:45 jm Exp $
-
- ********************************************************************/
-/**
- @file smallft.h
- @brief Discrete Rotational Fourier Transform (DRFT)
-*/
-
-#ifndef _V_SMFT_H_
-#define _V_SMFT_H_
-
-
-#ifdef __cplusplus
-extern "C" {-#endif
-
-/** Discrete Rotational Fourier Transform lookup */
-struct drft_lookup{- int n;
- float *trigcache;
- int *splitcache;
-};
-
-extern void spx_drft_forward(struct drft_lookup *l,float *data);
-extern void spx_drft_backward(struct drft_lookup *l,float *data);
-extern void spx_drft_init(struct drft_lookup *l,int n);
-extern void spx_drft_clear(struct drft_lookup *l);
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif
--- a/libcelt/testcelt.c
+++ b/libcelt/testcelt.c
@@ -97,6 +97,8 @@
celt_int16_t in[frame_size*channels];
celt_int16_t out[frame_size*channels];
fread(in, sizeof(short), frame_size*channels, fin);
+ if (feof(fin))
+ break;
len = celt_encode(enc, in, data, bytes_per_packet);
if (len <= 0)
{