function: PCM data vector blocking, windowing and dis/reassembly
author: Monty <xiphmont@mit.edu>
modifications by: Monty
- last modification date: Jul 28 1999
+ last modification date: Jul 29 1999
Handle windowing, overlap-add, etc of the PCM vectors. This is made
more amusing by Vorbis' current two allowed block sizes (512 and 2048
#include "envelope.h"
#include "mdct.h"
-/* pcm accumulator and multipliers
- examples (not exhaustive):
+/* pcm accumulator examples (not exhaustive):
- <-------------- lW----------------->
+ <-------------- lW ---------------->
<--------------- W ---------------->
: .....|..... _______________ |
: .''' | '''_--- | |\ |
|<------ Sl ------>| > Sr < |endW
|beginSl |endSl | |endSr
|beginW |endlW |beginSr
- mult[0] mult[n]
|< lW >|
|beginW | |endlW
mult[0] |beginSl mult[n]
- <-------------- lW----------------->
- |<-W-->|
-: .............. __ | |
-: .''' |`/ \ | |
-:.....''' |/`...\|...|
-:.........................|__||__|___|
- |Sl||Sr|endW
- | || |endSr
- | ||beginSr
- | |endSl
+ <-------------- lW ----------------->
+ |<--W-->|
+: .............. ___ | |
+: .''' |`/ \ | |
+:.....''' |/`....\|...|
+:.........................|___|___|___|
+ |Sl |Sr |endW
+ | | |endSr
+ | |beginSr
+ | |endSl
|beginSl
|beginW
- mult[0]
- mult[n]
*/
static int _vds_shared_init(vorbis_dsp_state *v,vorbis_info *vi){
--- /dev/null
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
+ * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
+ * PLEASE READ THESE TERMS DISTRIBUTING. *
+ * *
+ * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999 *
+ * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company *
+ * http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: LPC low level routines
+ author: Monty <monty@xiph.org>
+ modifications by: Monty
+ last modification date: Aug 02 1999
+
+ ********************************************************************/
+
+/* Some of these routines (autocorrelator, LPC coefficient estimator)
+ are directly derived from and/or modified from code written by
+ Jutta Degener and Carsten Bormann; thus we include their copyright
+ below. The entirety of this file is freely redistributable on the
+ condition that both of these copyright notices are preserved
+ without modification. */
+
+/* Preserved Copyright: *********************************************/
+
+/* Copyright 1992, 1993, 1994 by Jutta Degener and Carsten Bormann,
+Technische Universita"t Berlin
+
+Any use of this software is permitted provided that this notice is not
+removed and that neither the authors nor the Technische Universita"t
+Berlin are deemed to have made any representations as to the
+suitability of this software for any purpose nor are held responsible
+for any defects of this software. THERE IS ABSOLUTELY NO WARRANTY FOR
+THIS SOFTWARE.
+
+As a matter of courtesy, the authors request to be informed about uses
+this software has found, about bugs in this software, and about any
+improvements that may be of general interest.
+
+Berlin, 28.11.1994
+Jutta Degener
+Carsten Bormann
+
+*********************************************************************/
+
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include "smallft.h"
+#include "lpc.h"
+
+/* This is pared down for Vorbis, where we only use LPC to encode
+ spectral envelope curves. Thus we only are interested in
+ generating the coefficients and recovering the curve from the
+ coefficients. Autocorrelation LPC coeff generation algorithm
+ invented by N. Levinson in 1947, modified by J. Durbin in 1959. */
+
+/* Input : n element envelope curve
+ Output: m lpc coefficients, excitation energy */
+
+double vorbis_curve_to_lpc(double *curve,int n,double *lpc,int m){
+ double aut[m+1],work[n+n],error;
+ int i,j;
+
+ /* input is a real curve. make it complex-real */
+ for(i=0;i<n;i++){
+ work[i*2]=curve[i];
+ work[i*2+1]=0;
+ }
+
+ n*=2;
+ fft_backward(n,work,NULL,NULL);
+
+ /* The autocorrelation will not be circular. Shift, else we lose
+ most of the power in the edges. */
+
+ for(i=0,j=n/2;i<n/2;){
+ double temp=work[i];
+ work[i++]=work[j];
+ work[j++]=temp;
+ }
+
+ /* autocorrelation, p+1 lag coefficients */
+
+ j=m+1;
+ while(j--){
+ double d=0;
+ for(i=j;i<n;i++)d+=work[i]*work[i-j];
+ aut[j]=d;
+ }
+
+ /* Generate lpc coefficients from autocorr values */
+
+ error=aut[0];
+ if(error==0){
+ memset(lpc,0,m*sizeof(double));
+ return 0;
+ }
+
+ for(i=0;i<m;i++){
+ double r=-aut[i+1];
+
+ /* Sum up this iteration's reflection coefficient; note that in
+ Vorbis we don't save it. If anyone wants to recycle this code
+ and needs reflection coefficients, save the results of 'r' from
+ each iteration. */
+
+ for(j=0;j<i;j++)r-=lpc[j]*aut[i-j];
+ r/=error;
+
+ /* Update LPC coefficients and total error */
+
+ lpc[i]=r;
+ for(j=0;j<i/2;j++){
+ double tmp=lpc[j];
+ lpc[j]+=r*lpc[i-1-j];
+ lpc[i-1-j]+=r*tmp;
+ }
+ if(i%2)lpc[j]+=lpc[j]*r;
+
+ error*=1.0-r*r;
+ }
+
+ /* we need the error value to know how big an impulse to hit the
+ filter with later */
+
+ return error;
+}
+
+/* One can do this the long way by generating the transfer function in
+ the time domain and taking the forward FFT of the result. The
+ results from direct calculation are cleaner and faster, however */
+
+static double vorbis_lpc_magnitude(double w,double *lpc, int m){
+ int k;
+ double real=1,imag=0;
+ for(k=0;k<m;k++){
+ real+=lpc[k]*cos((k+1)*w);
+ imag+=lpc[k]*sin((k+1)*w);
+ }
+ return(1/sqrt(real*real+imag*imag));
+}
+
+/* generate the whole freq response curve on an LPC IIR filter */
+
+void vorbis_lpc_to_curve(double *curve,int n,double *lpc, double amp,int m){
+ int i;
+ double w=1./n*M_PI;
+ for(i=0;i<n;i++)
+ curve[i]=vorbis_lpc_magnitude(i*w,lpc,m)*amp;
+}
+
+/* find frequency response of LPC filter only at nonsero residue
+ points and apply the envelope to the residue */
+
+void vorbis_lpc_apply(double *residue,int n,double *lpc, double amp,int m){
+ int i;
+ double w=1./n*M_PI;
+ for(i=0;i<n;i++)
+ if(residue[i])
+ residue[i]*=vorbis_lpc_magnitude(i*w,lpc,m)*amp;
+}
+
+
--- /dev/null
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
+ * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
+ * PLEASE READ THESE TERMS DISTRIBUTING. *
+ * *
+ * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999 *
+ * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company *
+ * http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: LPC low level routines
+
+ ********************************************************************/
+
+#ifndef _V_LPC_H_
+#define _V_LPC_H_
+
+extern double vorbis_curve_to_lpc(double *curve,int n,double *lpc,int m);
+extern void vorbis_lpc_to_curve(double *curve,int n,double *lpc,
+ double amp,int m);
+extern void vorbis_lpc_apply(double *residue,int n,double *lpc,
+ double amp,int m);
+
+#endif
--- /dev/null
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
+ * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
+ * PLEASE READ THESE TERMS DISTRIBUTING. *
+ * *
+ * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999 *
+ * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company *
+ * http://www.xiph.org/ *
+ * *
+ ********************************************************************
+
+ function: LSP (also called LSF) conversion routines
+ author: Monty <monty@xiph.org>
+ modifications by: Monty
+ last modification date: Aug 03 1999
+
+ The LSP generation code is taken (with minimal modification) from
+ "On the Computation of the LSP Frequencies" by Joseph Rothweiler
+ <rothwlr@altavista.net>, available at:
+
+ http://www2.xtdl.com/~rothwlr/lsfpaper/lsfpage.html
+
+ ********************************************************************/
+
+#include <math.h>
+#include <string.h>
+#include <stdlib.h>
+
+void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m){
+ int i,j,m2=m/2;
+ double O[m2],E[m2];
+ double A,Ae[m2+1],Ao[m2+1];
+ double B,Be[m2],Bo[m2];
+ double temp;
+
+ /* even/odd roots setup */
+ for(i=0;i<m2;i++){
+ O[i]=-2.*cos(lsp[i*2]);
+ E[i]=-2.*cos(lsp[i*2+1]);
+ }
+
+ /* set up impulse response */
+ for(j=0;j<m2;j++){
+ Ae[j]=0.;
+ Ao[j]=1.;
+ Be[j]=0.;
+ Bo[j]=1.;
+ }
+ Ao[j]=1.;
+ Ae[j]=1.;
+
+ /* run impulse response */
+ for(i=1;i<m+1;i++){
+ A=B=0.;
+ for(j=0;j<m2;j++){
+ temp=O[j]*Ao[j]+Ae[j];
+ Ae[j]=Ao[j];
+ Ao[j]=A;
+ A+=temp;
+
+ temp=E[j]*Bo[j]+Be[j];
+ Be[j]=Bo[j];
+ Bo[j]=B;
+ B+=temp;
+ }
+ lpc[i-1]=(A+Ao[j]+B-Ae[j])/2;
+ Ao[j]=A;
+ Ae[j]=B;
+ }
+}
+
+static void kw(double *r,int n) {
+ double s[n/2+1];
+ double c[n+1];
+ int i, j, k;
+
+ s[0] = 1.0;
+ s[1] = -2.0;
+ s[2] = 2.0;
+ for(i=3;i<=n/2;i++) s[i] = s[i-2];
+
+ for(k=0;k<=n;k++) {
+ c[k] = r[k];
+ j = 1;
+ for(i=k+2;i<=n;i+=2) {
+ c[k] += s[j]*r[i];
+ s[j] -= s[j-1];
+ j++;
+ }
+ }
+ for(k=0;k<=n;k++) r[k] = c[k];
+}
+
+
+static int comp(const void *a,const void *b){
+ return(*(double *)a<*(double *)b);
+}
+
+/* CACM algorithm 283. */
+static void cacm283(double *a,int ord,double *r){
+ int i, k;
+ double val, p, delta, error;
+ double rooti;
+
+ for(i=0; i<ord;i++) r[i] = 2.0 * (i+0.5) / ord - 1.0;
+
+ for(error=1 ; error > 1.e-12; ) {
+ error = 0;
+ for( i=0; i<ord; i++) { /* Update each point. */
+ rooti = r[i];
+ val = a[ord];
+ p = a[ord];
+ for(k=ord-1; k>= 0; k--) {
+ val = val * rooti + a[k];
+ if (k != i) p *= rooti - r[k];
+ }
+ delta = val/p;
+ r[i] -= delta;
+ error += delta*delta;
+ }
+ }
+
+ /* Replaced the original bubble sort with a real sort. With your
+ help, we can eliminate the bubble sort in our lifetime. --Monty */
+
+ qsort(r,ord,sizeof(double),comp);
+
+}
+
+/* Convert lpc coefficients to lsp coefficients */
+void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m){
+ int order2=m/2;
+ double g1[order2+1], g2[order2+1];
+ double g1r[order2+1], g2r[order2+1];
+ int i;
+
+ /* Compute the lengths of the x polynomials. */
+ /* Compute the first half of K & R F1 & F2 polynomials. */
+ /* Compute half of the symmetric and antisymmetric polynomials. */
+ /* Remove the roots at +1 and -1. */
+
+ g1[order2] = 1.0;
+ for(i=0;i<order2;i++) g1[order2-i-1] = lpc[i]+lpc[m-i-1];
+ g2[order2] = 1.0;
+ for(i=0;i<order2;i++) g2[order2-i-1] = lpc[i]-lpc[m-i-1];
+
+ for(i=0; i<order2;i++) g1[order2-i-1] -= g1[order2-i];
+ for(i=0; i<order2;i++) g2[order2-i-1] += g2[order2-i];
+
+ /* Convert into polynomials in cos(alpha) */
+ kw(g1,order2);
+ kw(g2,order2);
+
+ /* Find the roots of the 2 even polynomials.*/
+
+ cacm283(g1,order2,g1r);
+ cacm283(g2,order2,g2r);
+
+ for(i=0;i<m;i+=2){
+ lsp[i] = acos(g1r[i/2]*.5);
+ lsp[i+1] = acos(g2r[i/2]*.5);
+ }
+}
+
--- /dev/null
+/********************************************************************
+ * *
+ * THIS FILE IS PART OF THE Ogg Vorbis SOFTWARE CODEC SOURCE CODE. *
+ * USE, DISTRIBUTION AND REPRODUCTION OF THIS SOURCE IS GOVERNED BY *
+ * THE GNU PUBLIC LICENSE 2, WHICH IS INCLUDED WITH THIS SOURCE. *
+ * PLEASE READ THESE TERMS DISTRIBUTING. *
+ * *
+ * THE OggSQUISH SOURCE CODE IS (C) COPYRIGHT 1994-1999 *
+ * by 1999 Monty <monty@xiph.org> and The XIPHOPHORUS Company *
+ * http://www.xiph.org/ *
+ * *
+ ********************************************************************/
+
+#ifndef _V_LSP_H_
+#define _V_LSP_H_
+
+extern void vorbis_lsp_to_lpc(double *lsp,double *lpc,int m);
+extern void vorbis_lpc_to_lsp(double *lpc,double *lsp,int m);
+
+#endif
--- /dev/null
+/******************************************************************
+ * CopyPolicy: GNU Public License 2 applies
+ * Copyright (C) 1998 Monty xiphmont@mit.edu, monty@xiph.org
+ *
+ * FFT implementation from OggSquish, minus cosine transforms,
+ * minus all but radix 2/4 case. In Vorbis we only need this
+ * cut-down version (used by the encoder, not decoder).
+ *
+ * To do more than just power-of-two sized vectors, see the full
+ * version I wrote for NetLib.
+ *
+ ******************************************************************/
+
+#include <stdlib.h>
+#include <math.h>
+#include "smallft.h"
+
+static void drfti1(int n, double *wa, int *ifac){
+ static int ntryh[4] = { 4,2,3,5 };
+ static double tpi = 6.28318530717958647692528676655900577;
+ double 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=(double)ld*argh;
+ fi=0.;
+ for (ii=2;ii<ido;ii+=2){
+ fi+=1.;
+ arg=fi*argld;
+ wa[i++]=cos(arg);
+ wa[i++]=sin(arg);
+ }
+ is+=ido;
+ }
+ l1=l2;
+ }
+}
+
+static void fdrffti(int n, double *wsave, int *ifac){
+
+ if (n == 1) return;
+ drfti1(n, wsave+n, ifac);
+}
+
+static void dradf2(int ido,int l1,double *cc,double *ch,double *wa1){
+ int i,k;
+ double 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,double *cc,double *ch,double *wa1,
+ double *wa2,double *wa3){
+ static double hsqt2 = .70710678118654752440084436210485;
+ int i,k,t0,t1,t2,t3,t4,t5,t6;
+ double 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 drftf1(int n,double *c,double *ch,double *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:
+ return; /* We're restricted to powers of two. just fail */
+
+ L110:
+ l2=l1;
+ }
+
+ if(na==1)return;
+
+ for(i=0;i<n;i++)c[i]=ch[i];
+}
+
+static void fdrfftf(int n,double *r,double *wsave,int *ifac){
+ if(n==1)return;
+ drftf1(n,r,wsave,wsave+n,ifac);
+}
+
+static void dradb2(int ido,int l1,double *cc,double *ch,double *wa1){
+ int i,k,t0,t1,t2,t3,t4,t5,t6;
+ double 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 dradb4(int ido,int l1,double *cc,double *ch,double *wa1,
+ double *wa2,double *wa3){
+ static double sqrt2=1.4142135623730950488016887242097;
+ int i,k,t0,t1,t2,t3,t4,t5,t6,t7,t8;
+ double 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 drftb1(int n, double *c, double *ch, double *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:
+ return; /* silently fail. we only do powers of two in this version */
+
+ L115:
+ l1=l2;
+ iw+=(ip-1)*ido;
+ }
+
+ if(na==0)return;
+
+ for(i=0;i<n;i++)c[i]=ch[i];
+}
+
+static void fdrfftb(int n, double *r, double *wsave, int *ifac){
+ if (n == 1)return;
+ drftb1(n, r, wsave, wsave+n, ifac);
+}
+
+void fft_forward(int n, double *buf,double *trigcache,int *splitcache){
+ int flag=0;
+
+ if(!trigcache || !splitcache){
+ trigcache=calloc(3*n,sizeof(double));
+ splitcache=calloc(32,sizeof(int));
+ fdrffti(n, trigcache, splitcache);
+ flag=1;
+ }
+
+ fdrfftf(n, buf, trigcache, splitcache);
+
+ if(flag){
+ free(trigcache);
+ free(splitcache);
+ }
+}
+
+void fft_backward(int n, double *buf, double *trigcache,int *splitcache){
+ int i;
+ int flag=0;
+
+ if(!trigcache || !splitcache){
+ trigcache=calloc(3*n,sizeof(double));
+ splitcache=calloc(32,sizeof(int));
+ fdrffti(n, trigcache, splitcache);
+ flag=1;
+ }
+
+ fdrfftb(n, buf, trigcache, splitcache);
+
+ for(i=0;i<n;i++)buf[i]/=n;
+
+ if(flag){
+ free(trigcache);
+ free(splitcache);
+ }
+}
+
+void fft_i(int n, double **trigcache, int **splitcache){
+ *trigcache=calloc(3*n,sizeof(double));
+ *splitcache=calloc(32,sizeof(int));
+ fdrffti(n, *trigcache, *splitcache);
+}
--- /dev/null
+/******************************************************************
+ * CopyPolicy: GNU Public License 2 applies
+ * Copyright (C) 1998 Monty xiphmont@mit.edu, monty@xiph.org
+ *
+ * FFT implementation from OggSquish, minus cosine transforms.
+ * Only convenience functions exposed
+ *
+ ******************************************************************/
+
+#ifndef _V_SMFT_H_
+#define _V_SMFT_H_
+
+extern void fft_forward(int n, double *buf, double *trigcache, int *sp);
+extern void fft_backward(int n, double *buf, double *trigcache, int *sp);
+extern void fft_i(int n, double **trigcache, int **splitcache);
+
+#endif
projects / platform / upstream / libvorbis.git / commitdiff