0.9 Appendix 2 - ffts with precomputed luts

Computing the fast fourier Page 1 / 1

This Appendix contains source code listings corresponding to the FFT implementations with precomputed coefficients in Implementation details .

#include <math.h>#include <complex.h>#include <stdio.h>#include <stdlib.h> 
typedef complex float data_t; 
#define W(N,k) (cexp(-2.0f * M_PI * I * (float)k / (float)N)) 
data_t *LUT; 
void ditfft2(data_t *in, data_t *out, int log2stride, int stride, int N) {if(N == 2) {
		out[0]   = in[0] + in[stride];
		out[N/2] = in[0] - in[stride];
	}else{ditfft2(in, out, log2stride+1, stride << 1, N >> 1);
		ditfft2(in+stride, out+N/2, log2stride+1, stride << 1, N >> 1);{	 /* k=0 -> no multiplication */
			data_t Ek = out[0];
			data_t Ok = out[N/2];
			out[0]   = Ek + Ok;
			out[N/2] = Ek - Ok;
		} 
		int k;for(k=1;k<N/2;k++) {
			data_t Ek = out[k];
			data_t Ok = out[(k+N/2)];
			data_t w = LUT[k<<log2stride];out[k]        = Ek + w * Ok;out[(k+N/2) ] = Ek - w * Ok;}
	}}
 void fft_init(int N) {
	LUT = malloc(N/2 * sizeof(data_t));int i;
	for(i=0;i<N/2;i++) LUT[i] = W(N,i);}
Simple radix-2 FFT with precomputed LUT

#include <complex.h>#include <stdio.h>#include <stdlib.h> 
typedef complex float data_t; 
#define W(N,k) (cexp(-2.0f * M_PI * I * (float)k / (float)N))data_t *LUT1;
data_t *LUT3; 
void splitfft(data_t *in, data_t *out, int log2stride, int stride, int N) {if(N == 1) {
		out[0] = in[0];}else if(N == 2) {
		out[0]   = in[0] + in[stride];
		out[N/2] = in[0] - in[stride];
	}else{splitfft(in, out, log2stride+1, stride << 1, N >> 1);
		splitfft(in+stride, out+N/2, log2stride+2, stride << 2, N >> 2);
	  splitfft(in+3*stride, out+3*N/4, log2stride+2, stride << 2, N >> 2);{	
			data_t Uk  = out[0];
			data_t Zk  = out[0+N/2];
			data_t Uk2 = out[0+N/4];
			data_t Zdk = out[0+3*N/4];
			out[0]       = Uk  + (Zk + Zdk);
			out[0+N/2]   = Uk  - (Zk + Zdk);
			out[0+N/4]   = Uk2 - I*(Zk - Zdk);
			out[0+3*N/4] = Uk2 + I*(Zk - Zdk);
		}int k;
		for(k=1;k<N/4;k++) {
			data_t Uk  = out[k];
			data_t Zk  = out[k+N/2];
			data_t Uk2 = out[k+N/4];
			data_t Zdk = out[k+3*N/4];
			data_t w1 = LUT1[k<<log2stride];data_t w3 = LUT3[k<<log2stride];out[k]       = Uk  + (w1*Zk + w3*Zdk);out[k+N/2]   = Uk  - (w1*Zk + w3*Zdk);out[k+N/4]   = Uk2 - I*(w1*Zk - w3*Zdk);out[k+3*N/4] = Uk2 + I*(w1*Zk - w3*Zdk);}
	}}
 void fft_init(int N) {
	LUT1 = malloc(N/4 * sizeof(data_t));LUT3 = malloc(N/4 * sizeof(data_t));
	int i;for(i=0;i<N/4;i++) LUT1[i] = W(N,i);for(i=0;i<N/4;i++) LUT3[i] = W(N,3*i);}
Simple split-radix FFT with precomputed LUT

#include <math.h>#include <complex.h>#include <stdio.h>#include <stdlib.h> 
typedef complex float data_t; 
#define W(N,k) (cexp(-2.0f * M_PI * I * (float)k / (float)N))data_t *LUT;
 void conjfft(data_t *base, int TN,
	data_t *in, data_t *out, int log2stride, int stride, int N) {if(N == 1) {
		if(in < base) in += TN;
		out[0] = in[0];}else if(N == 2) {
		data_t *i0 = in, *i1 = in + stride;if(i0 < base) i0 += TN;
		if(i1 < base) i1 += TN;
		out[0]   = *i0 + *i1;
		out[N/2] = *i0 - *i1;
	}else{conjfft(base, TN, in, out, log2stride+1, stride << 1, N >> 1);
		conjfft(base, TN, in+stride, out+N/2, log2stride+2, stride << 2, N >> 2);
	  conjfft(base, TN, in-stride, out+3*N/4, log2stride+2, stride << 2, N >> 2);{	
			data_t Uk  = out[0];
			data_t Zk  = out[0+N/2];
			data_t Uk2 = out[0+N/4];
			data_t Zdk = out[0+3*N/4];
			out[0]       = Uk  + (Zk + Zdk);
			out[0+N/2]   = Uk  - (Zk + Zdk);
			out[0+N/4]   = Uk2 - I*(Zk - Zdk);
			out[0+3*N/4] = Uk2 + I*(Zk - Zdk);
		}int k;
		for(k=1;k<N/4;k++) {
			data_t Uk  = out[k];
			data_t Zk  = out[k+N/2];
			data_t Uk2 = out[k+N/4];
			data_t Zdk = out[k+3*N/4];
			data_t w = LUT[k<<log2stride];out[k]       = Uk  + (w*Zk + conj(w)*Zdk);out[k+N/2]   = Uk  - (w*Zk + conj(w)*Zdk);out[k+N/4]   = Uk2 - I*(w*Zk - conj(w)*Zdk);out[k+3*N/4] = Uk2 + I*(w*Zk - conj(w)*Zdk);}
	}}void fft_init(int N) {
	LUT = malloc(N/4 * sizeof(data_t));int i;
	for(i=0;i<N/4;i++) LUT[i] = W(N,i);}
Simple conjugate-pair FFT with precomputed LUT

#include <math.h>#include <complex.h>#include <stdio.h>#include <stdlib.h> 
typedef complex float data_t; 
#define W(N,k) (cexp(-2.0f * M_PI * I * (float)(k) / (float)(N))) 
floats(int n, int k) {
  if (n <= 4) return 1.0f;
   int k4 = k % (n/4);
   if (k4 <= n/8)
		return (s(n/4,k4) * cosf(2.0f * M_PI * (float)k4 / (float)n));  return (s(n/4,k4) * sinf(2.0f * M_PI * (float)k4 / (float)n));
} 
data_t *LUT, *LUT0, *LUT1, *LUT2;float *s2, *s4;
 void tangentfft8(data_t *base, int TN, data_t *in, data_t *out, int log2stride, 	
		int stride, int N) {  if(N == 1) {
		if(in < base) in += TN;
		out[0] = in[0];  }else if(N == 2) {
		data_t *i0 = in, *i1 = in + stride;if(i0 < base) i0 += TN;
		if(i1 < base) i1 += TN;
		out[0]   = *i0 + *i1;
		out[N/2] = *i0 - *i1;
  }else if(N == 4) {    tangentfft8(base, TN, in, out, log2stride+1, stride << 1, N >> 1);
    tangentfft8(base, TN, in+stride, out+2, log2stride+1, stride << 1, N >> 1);
     data_t temp1 = out[0] + out[2];
    data_t temp2 = out[0] - out[2];    out[0] = temp1;    out[2] = temp2;    temp1 = out[1] - I*out[3];
    temp2 = out[1] + I*out[3];    out[1] = temp1;    out[3] = temp2; 
  }else{    tangentfft8(base, TN, in, out, log2stride+2, stride << 2, N >> 2);
    tangentfft8(base, TN, in+(stride*2), out+2*N/8, log2stride+3, stride << 3, N >> 3);
    tangentfft8(base, TN, in-(stride*2), out+3*N/8, log2stride+3, stride << 3, N >> 3);
    tangentfft8(base, TN, in+(stride), out+4*N/8, log2stride+2, stride << 2, N >> 2);
    tangentfft8(base, TN, in-(stride), out+6*N/8, log2stride+2, stride << 2, N >> 2);
		int k;    for(k=0;k<N/8;k++) {
			data_t w0 = LUT0[k<<log2stride];data_t w1 = LUT1[k<<log2stride];data_t w2 = LUT2[k<<log2stride]; 
      data_t zk_p   = w0       * out[k+4*N/8];
      data_t zk_n   = conj(w0) * out[k+6*N/8];
      data_t zk2_p  = w1       * out[k+5*N/8];
      data_t zk2_n  = conj(w1) * out[k+7*N/8];
      data_t uk     = out[k]                  * s4[k<<log2stride];      data_t uk2    = out[k+N/8]              * s4[k+N/8 << log2stride];      data_t yk_p   = w2       * out[k+2*N/8];      data_t yk_n   = conj(w2) * out[k+3*N/8];data_t y0 = (yk_p + yk_n)*s2[k<<log2stride];data_t y1 = (yk_p - yk_n)*I*s2[k+N/8 << log2stride]; 
      out[k]       = uk + y0 + (zk_p + zk_n);
      out[k+4*N/8] = uk + y0 - (zk_p + zk_n);
      out[k+2*N/8] = uk - y0 - I*(zk_p - zk_n);
      out[k+6*N/8] = uk - y0 + I*(zk_p - zk_n);
      out[k+1*N/8] = uk2 - y1 +   (zk2_p + zk2_n);
      out[k+3*N/8] = uk2 + y1 - I*(zk2_p - zk2_n);
      out[k+5*N/8] = uk2 - y1 -   (zk2_p + zk2_n);
      out[k+7*N/8] = uk2 + y1 + I*(zk2_p - zk2_n);
    }  }
 }
 void tangentfft4(data_t *base, int TN, data_t *in, data_t *out, int log2stride,
		int stride, int N) {if(N == 1) {
		if(in < base) in += TN;
		out[0] = in[0];}else if(N == 2) {
		data_t *i0 = in, *i1 = in + stride;if(i0 < base) i0 += TN;
		if(i1 < base) i1 += TN;
		out[0]   = *i0 + *i1;
		out[N/2] = *i0 - *i1;
	}else{tangentfft4(base, TN, in, out, log2stride+1, stride << 1, N >> 1);
		tangentfft8(base, TN, in+stride, out+N/2, log2stride+2, stride << 2, N >> 2);
	  tangentfft8(base, TN, in-stride, out+3*N/4, log2stride+2, stride << 2, N >> 2);{	
			data_t Uk  = out[0];
			data_t Zk  = out[0+N/2];
			data_t Uk2 = out[0+N/4];
			data_t Zdk = out[0+3*N/4];
			out[0]       = Uk  + (Zk + Zdk);
			out[0+N/2]   = Uk  - (Zk + Zdk);
			out[0+N/4]   = Uk2 - I*(Zk - Zdk);
			out[0+3*N/4] = Uk2 + I*(Zk - Zdk);
		}int k;
		for(k=1;k<N/4;k++) {
			data_t Uk  = out[k];
			data_t Zk  = out[k+N/2];
			data_t Uk2 = out[k+N/4];
			data_t Zdk = out[k+3*N/4];
			data_t w = LUT[k<<log2stride];out[k]       = Uk  + (w*Zk + conj(w)*Zdk);out[k+N/2]   = Uk  - (w*Zk + conj(w)*Zdk);out[k+N/4]   = Uk2 - I*(w*Zk - conj(w)*Zdk);out[k+3*N/4] = Uk2 + I*(w*Zk - conj(w)*Zdk);}
	}}
 void fft_init(int N) {
	LUT0 = malloc(N/8 * sizeof(data_t));LUT1 = malloc(N/8 * sizeof(data_t));
	LUT2 = malloc(N/8 * sizeof(data_t)); 	LUT = malloc(N/4 * sizeof(data_t));
 s2 = malloc(N/4 * sizeof(float));
	s4 = malloc(N/4 * sizeof(float)); 
	int i;for(i=0;i<N/8;i++) LUT0[i] = W(N,i)*s(N/4,i)/s(N,i);for(i=0;i<N/8;i++) LUT1[i] = W(N,i+N/8)*s(N/4,i+N/8)/s(N,i+N/8);for(i=0;i<N/8;i++) LUT2[i] = W(N,2*i)*s(N/8,i)/s(N/2,i);for(i=0;i<N/4;i++) LUT[i] = W(N,i)*s(N/4,i);for(i=0;i<N/4;i++) s4[i] = s(N/4,i)/s(N,i);for(i=0;i<N/4;i++) s2[i] = s(N/2,i)/s(N,i);}
Simple tangent FFT with precomputed LUT

Questions & Answers

A golfer on a fairway is 70 m away from the green, which sits below the level of the fairway by 20 m. If the golfer hits the ball at an angle of 40° with an initial speed of 20 m/s, how close to the green does she come?

Aislinn Reply

tijani

what is titration

John Reply

what is physics

Siyaka Reply

A mouse of mass 200 g falls 100 m down a vertical mine shaft and lands at the bottom with a speed of 8.0 m/s. During its fall, how much work is done on the mouse by air resistance

Jude Reply

Can you compute that for me. Ty

Jude

what is the dimension formula of energy?

David Reply

what is viscosity?

David

what is inorganic

emma Reply

what is chemistry

Youesf Reply

what is inorganic

emma

Chemistry is a branch of science that deals with the study of matter,it composition,it structure and the changes it undergoes

Adjei

please, I'm a physics student and I need help in physics

Adjanou

chemistry could also be understood like the sexual attraction/repulsion of the male and female elements. the reaction varies depending on the energy differences of each given gender. + masculine -female.

Pedro

A ball is thrown straight up.it passes a 2.0m high window 7.50 m off the ground on it path up and takes 1.30 s to go past the window.what was the ball initial velocity

Krampah Reply

2. A sled plus passenger with total mass 50 kg is pulled 20 m across the snow (0.20) at constant velocity by a force directed 25° above the horizontal. Calculate (a) the work of the applied force, (b) the work of friction, and (c) the total work.

Sahid Reply

you have been hired as an espert witness in a court case involving an automobile accident. the accident involved car A of mass 1500kg which crashed into stationary car B of mass 1100kg. the driver of car A applied his brakes 15 m before he skidded and crashed into car B. after the collision, car A s

Samuel Reply

can someone explain to me, an ignorant high school student, why the trend of the graph doesn't follow the fact that the higher frequency a sound wave is, the more power it is, hence, making me think the phons output would follow this general trend?

Joseph Reply

Nevermind i just realied that the graph is the phons output for a person with normal hearing and not just the phons output of the sound waves power, I should read the entire thing next time

Joseph

Follow up question, does anyone know where I can find a graph that accuretly depicts the actual relative "power" output of sound over its frequency instead of just humans hearing

Joseph

"Generation of electrical energy from sound energy | IEEE Conference Publication | IEEE Xplore" ***ieeexplore.ieee.org/document/7150687?reload=true

Ryan

what's motion

Maurice Reply

what are the types of wave

Maurice

answer

Magreth

progressive wave

Magreth

hello friend how are you

Muhammad Reply

fine, how about you?

Mohammed

Mujahid

A string is 3.00 m long with a mass of 5.00 g. The string is held taut with a tension of 500.00 N applied to the string. A pulse is sent down the string. How long does it take the pulse to travel the 3.00 m of the string?

yasuo Reply

Who can show me the full solution in this problem?

Reofrir Reply

Got questions? Join the online conversation and get instant answers!

Jobilize.com Reply

<< Chapter < Page Page > Chapter >>

Read also:

Get Jobilize Job Search Mobile App in your pocket Now!

100% Free Mobile Applications
Receive real-time job alerts and never miss the right job again

Source: OpenStax, Computing the fast fourier transform on simd microprocessors. OpenStax CNX. Jul 15, 2012 Download for free at http://cnx.org/content/col11438/1.2

Google Play and the Google Play logo are trademarks of Google Inc.

Notification Switch

Would you like to follow the 'Computing the fast fourier transform on simd microprocessors' conversation and receive update notifications?

Ask

	8 Psychology MCQ 2011 1 Exam By John Gabrieli Start Exam
	30 Biology 30 Plant Form and Physiology MCQ By OpenStax Start Quiz
©flickr: Donnie	How well do you know Beanie Boos? By Angela Eckman Start Quiz
	NCE Ch 04 Social and Cultural Foundations By Anh Dao Start Quiz
©flickr: Abraham	3 Biology 3 By Sarah Warren Start Test
	Principles of microeconomics for ap® courses By OpenStax Read Online Course
	6 Physiotherapy Flashcards Set 6 By Rhodes Start Flashcards
	14 AP 14 Brain Cranial Nerves MCQ By OpenStax Start Quiz
	22 AP Key Terms 22 The Respiratory System By OpenStax Start Key Terms
	2 Understanding Societies 2 By Jessica Collett Start Exam