 Notice. New forum software under development. It's going to miss a few functions and look a bit ugly for a while, but I'm working on it full time now as the old forum was too unstable. Couple days, all good. If you notice any issues, please contact me. |
Forum Index : Microcontroller and PC projects : MM2(+) Fast Fourier Transform
| Author | Message | ||||
| matherp Guru  Joined: 11/12/2012 Location: United KingdomPosts: 8584 |
Please find attached a CFunction that does an in-place Fast Fourier Transform (FFT) of any dataset of size 2^n elements. The data points are specified as complex numbers however, only the real component needs providing (the imaginary component should be set to zero). The format of the provided data must be as a 2*n floating point array specified as DIM amplitudes!(1,n-1). The FFT subroutine returns the array transformed into the frequency domain as complex numbers. The magnitude of a given frequency can be calculated by sqr(real^2 + imaginary^2) The phase of a given frequency can be calculated by atan2(imaginary,real) Functions to do these calculations are given in this post The CFunction is a port of the Liberty Basic code on rosettacode.org. This was chosen as it doesn't use recursive coding which can give problems in CFunctions. On the small test dataset the Cfunction take less than 1msec to complete whereas the Basic takes 36msec. OPTION EXPLICIT OPTION DEFAULT NONE dim i% DIM amplitudes!(1,7)=(1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0) 'Test data set from rosettacode.org FFT 8,amplitudes!() for i%=0 to 7 print "Real=",amplitudes!(0,i%)," Imaginary=",amplitudes!(1,i%) next i% end ' ' In place FFT for 2^n datasets ' Parameters are: ' Number of data points (must be power of 2) ' Two dimensional array of floats dimensioned (1,N) where 0,N has the real component and 1,N has the imaginary component ' CSUB FFT 'In place FFT for 2^n datasets 00000068 'cos 27BDFFE0 3C039D00 AFB20018 8C62009C 00809021 3C043FC9 AFBF001C AFB10014 AFB00010 24840FDB 8C70006C 0040F809 8C710060 00402021 0220F809 02402821 00402021 0200C821 8FBF001C 8FB20018 8FB10014 8FB00010 03200008 27BD0020 'Re 000420C0 00A42021 03E00008 8C820000 'Im 000420C0 00A42021 03E00008 8C820004 'setRe 000420C0 00C42021 03E00008 AC850000 'setIm 000420C0 00C42021 03E00008 AC850004 'pflt 27BDFFC0 24020020 AFB00034 AFA20010 3C109D00 00801821 8E020034 AFB10038 27A40018 00A08821 AFBF003C 00602821 24060008 0040F809 24070004 8E02002C 0040F809 27A40018 16200008 8E020008 0040F809 2404000A 8FBF003C 8FB10038 8FB00034 03E00008 27BD0040 0040F809 24040020 8FBF003C 8FB10038 8FB00034 03E00008 27BD0040 'pint 27BDFFC8 2402000A AFB00030 AFA20010 3C109D00 00803021 8E020030 AFBF0034 AFA50028 27A40018 0040F809 00063FC3 8E02002C 0040F809 27A40018 8FA50028 14A00007 8E020008 0040F809 2404000A 8FBF0034 8FB00030 03E00008 27BD0038 0040F809 24040020 8FBF0034 8FB00030 03E00008 27BD0038 'main 27BDFF88 AFB40060 AFB3005C 8C940000 3C139D00 8E620040 AFBF0074 001420C0 AFBE0070 AFB7006C AFB60068 AFB50064 AFB20058 AFB10054 AFB00050 0040F809 00A08021 00408821 00141880 02231821 8E62009C 3C044049 AFA30034 0040F809 24840FDB 00409021 8E620080 00142FC3 02802021 8E770070 8E760064 8E74007C 0040F809 8E75005C 02E0F809 00402021 0040B821 8E62009C 3C044000 0040F809 8E7E0070 03C0F809 00402021 00402821 02C0F809 02E02021 0040B021 3C043EFF 8E62009C 0040F809 3484F2E5 00402821 02A0F809 02C02021 0280F809 00402021 2444FFFF 1880017D AFA20040 00001021 24130002 24420001 0044182A 1460FFFD 00139840 8FA20040 00131843 000227C2 00822021 00042043 00446823 00133083 25A5FFFF AFA30038 18A00170 AFA60018 24020002 AFA20010 00001021 8FA60010 24420001 00063040 0045182A 1460FFFB AFA60010 2484FFFF 18800162 00001021 24160002 24420001 0044182A 1460FFFD 0016B040 19A0001F AE200000 8FA50010 24090001 240C0001 00055FC2 01655821 8FA30010 000B5843 006B5023 014B102A 1440000E 000B1080 01652021 00053080 02201821 02221021 8C680000 00852021 00853823 01284021 0147382A AC480000 00661821 10E0FFF8 00461021 258C0001 01AC102A 14400003 00094840 1000FFE6 01602821 3C149D00 8E82009C 3C044000 8E950064 0040F809 8E970058 00402821 02E0F809 02402021 00409021 8E820080 02602021 0040F809 00132FC3 00402821 02A0F809 02402021 8FA40018 0480002B 0040A821 00009021 3C139D00 8E620080 02402021 00122FC3 0040F809 8E740058 00402021 0280F809 02A02821 00402021 0411FEFB 00000000 8FA50034 00121880 0040A021 00A31821 8E620060 00002021 02802821 0040F809 AC740000 8FA60038 8FA30034 00D22023 00042080 00642021 8E630060 02802821 AC820000 0060F809 00002021 8FA50038 8FA60034 02452021 8FA50018 00042080 26520001 00C42021 00B2182A 1060FFD9 AC820000 0220F021 0000B821 AFA00014 0220A821 10000004 00009021 0256102A 10400034 26B50004 8EB40000 8FC20000 72D41802 02F29821 0062A021 0274102A 1040FFF6 26520001 02802021 02002821 0411FEE1 00000000 02602021 02002821 AFA20048 0411FEDC 00000000 00402821 02802021 02003021 0411FEDF 00000000 8FA30048 02003021 00602821 02602021 0411FED9 00000000 02802021 02002821 0411FED1 00000000 02602021 02002821 AFA20048 0411FECC 00000000 00402821 02802021 02003021 0411FECF 00000000 8FA30048 02602021 00602821 02003021 0411FEC9 00000000 0256102A 1440FFCE 26B50004 8FA40014 8FA50010 24840001 0085102A AFA40014 02F6B821 1440FFC0 27DE0004 8FA60040 18C000AA 8FA30018 8FA50034 00031080 24040001 00A21021 AFA0003C AFA40020 AFA20044 3C119D00 8FA30038 8FA40020 8FA60034 0064001A 008001F4 8FA20044 00042840 AFA50028 AFA60018 AFA2001C AFA0002C 00001812 AFA30024 00031880 AFA30030 8FA40024 1880007B 8FA3002C 8FB2002C 0000A021 02402021 02002821 0411FE8F 00000000 8FA60020 02402021 02002821 02469821 AFA20010 0411FE8C 00000000 0040B821 8FA20018 02602021 8C560000 02002821 8E350058 0411FE80 00000000 00402821 02A0F809 02C02021 8FA3001C 02602021 8C760000 02002821 0040F021 8E350058 0411FE79 00000000 00402821 02A0F809 02C02021 8FA4001C 02002821 8C960000 02602021 8E350058 AFA20014 0411FE6A 00000000 00402821 02A0F809 02C02021 8FA50018 02602021 8CA30000 02002821 8E350058 0040B021 AFA30048 0411FE62 00000000 8FA30048 00402821 02A0F809 00602021 8E230060 0040A821 8FA40010 8E22005C 03C02821 0040F809 AFA30048 8FA30048 8FA50014 0060F809 00402021 00402821 02003021 02402021 0411FE52 00000000 8E23005C 02E02021 02C02821 0060F809 AFA30048 8FA30048 00402021 0060F809 02A02821 00402821 02003021 02402021 0411FE48 00000000 8E230060 8E22005C 8FA40010 8FA50014 0040F809 AFA30048 8FA30048 03C02821 0060F809 00402021 00402821 02003021 02602021 0411FE35 00000000 8E3E0060 02E02021 03C0F809 02C02821 00402021 03C0F809 02A02821 00402821 02003021 02602021 0411FE2D 00000000 8FA60028 8FA20024 26940001 1682FF8A 02469021 8FA3002C 8FA40020 24630001 8FA50018 8FA60030 AFA3002C 0064102A 8FA3001C 00A62821 00661821 AFA50018 1440FF78 AFA3001C 8FA4003C 8FA50040 24840001 0085102A 10400004 AFA4003C 8FA60028 1000FF60 AFA60020 8FBF0074 8FBE0070 8FB7006C 8FB60068 8FB50064 8FB40060 8FB3005C 8FB20058 8FB10054 8FB00050 03E00008 27BD0078 1000FE89 24130002 1000FEA3 24160002 24020002 1000FE98 AFA20010 End CSUB Basic equivalent: option explicit option default none DIM INTEGER S,P,R1,R2,R3,R4,R,i,j,S1,S2,P1,P2,DV,DP,HA,PT,IP,H,G,T,stage,D,Z,L,LS,A,B,N,X,M; dim float K,COX,temp,F1,F2,X1,X2,X3,X4,dummy,aa! P =8 Dim float Re( P), Im( P), Co( P) for N =0 to P -1 read dummy: Re( N) =dummy read dummy: Im( N) =dummy next N S =int( log( P) /log( 2) +0.9999) R1 =2^S R =R1 -1 R2 =div( R1, 2) R4 =div( R1, 4) R3 =R4 +R2 data 1, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0 S2 =div( S, 2) S1 =S -S2 P1 =2^S1 P2 =2^S2 dim integer V( P1 -1) V( 0) =0 DV =1 DP =P1 for J =1 to S1 HA =div( DP, 2) PT =P1 -HA for I =HA to PT step DP V( I) =V( I -HA) +DV next I DV =DV +DV DP =HA next J K =2 *Pi /R1 for X =0 to R4 COX =cos( K *X) Co( X) =COX Co( R2 -X) =0 -COX Co( R2 +X) =0 -COX next X ' print "FFT: bit reversal" for I =0 to P1 -1 IP =I *P2 for J =0 to P2 -1 H =IP +J G =V( J) *P2 +V( I) if G >H then temp =Re( G): Re( G) =Re( H): Re( H) =temp if G >H then temp =Im( G): Im( G) =Im( H): Im( H) =temp next J next I T =1 for stage =0 to S -1 ' print " Stage:- "; stage D =div( R2, T) for Z =0 to T -1 L =D *Z LS =L +R4 for I =0 to D -1 A =2 *I *T +Z B =A +T F1 =Re( A) F2 =Im( A) X1 =Co( L) *Re( B) X2 =Co( LS) *Im( B) X3 =Co( LS) *Re( B) X4 =Co( L) *Im( B) Re( A) =F1 +X1 -X2 Im( A) =F2 +X3 +X4 Re( B) =F1 -X1 +X2 Im( B) =F2 -X3 -X4 next I next Z T =T +T next stage print " M Re( M) Im( M)" for M =0 to R if abs( Re( M)) <10^-5 then Re( M) =0 if abs( Im( M)) <10^-5 then Im( M) =0 print " "; M, Re( M), Im( M) next M end wait function div( a as integer, b as integer) as integer div =a \b end function end C Source #define Version 100 //Version 1.00 #include "../cfunctions.h" float cos(float x){ return FSin(FSub(LoadFloat(0x3fc90fdb),x)); } float Re(int N, float amplitude[]){ return amplitude[N*2]; } float Im(int N, float amplitude[]){ return amplitude[N*2+1]; } void setRe(int N, float val, float amplitude[]){ amplitude[N*2]=val; } void setIm(int N, float val, float amplitude[]){ amplitude[N*2+1]=val; } void pflt(float p,int f){ char buff[20]; FloatToStr(buff,p,8,4,32); MMPrintString(buff); if(f){ putConsole(32); } else { putConsole(10); } } void pint(int p,int f){ char x[10]; IntToStr(x,p,10); MMPrintString(x); if(f){ putConsole(32); } else { putConsole(10); } } void main(long *n, float *amplitude){ int P=*n,R1,R2,R3,R4,R,l,j,S1,S2,P1,P2,DV,DP,HA,PT,IP,H,G,T,stage,D,Z,L,LS,A,B; int *V=GetTempMemory(P * 8); float *Co; Co=(void*)V+P*4;; float pi=LoadFloat(0x40490FDB),K,COX,temp,F1,F2,X1,X2,X3,X4; int S=FloatToInt(FAdd(FDiv(FLog(IntToFloat(P)),FLog(LoadFloat(0x40000000))),LoadFloat(0x3EFFF2E5))); R1=2; for(l=0;l<S-1;l++)R1*=2; R=R1-1; R2=R1/2; R4=R1/4; R3=R4+R2; S2=S/2; S1=S-S2; P1=2; for(l=0;l<S1-1;l++)P1*=2; P2=2; for(l=0;l<S2-1;l++)P2*=2; V[0]=0; DV=1; DP=P1; for(j=1;j<=S1;j++){ HA=DP/2; PT=P1-HA; for(l=HA;l<=PT;l+=DP)V  =V[l-HA]+DV;
DV*=2; DP=HA; } K=FDiv(FMul(pi,LoadFloat(0x40000000)),IntToFloat(R1)); for(l=0;l<=R4;l++){ COX=cos(FMul(IntToFloat(l),K)); Co =COX;
Co[R2-l]=FSub(0,COX); Co[R2+l]=FSub(0,COX); } for(l=0;l<P1;l++){ IP=l*P2; for(j=0;j<P2;j++){ H=IP+j; G=V[j]*P2+V ;
if(G>H){ temp=Re(G,amplitude); setRe(G,Re(H,amplitude),amplitude); setRe(H,temp,amplitude); temp=Im(G,amplitude); setIm(G,Im(H,amplitude),amplitude); setIm(H,temp,amplitude); } } } T=1; for(stage=0;stage<S;stage++){ D=R2/T; for(Z=0;Z<T;Z++){ L=D*Z; LS=L+R4; for(l=0;l<D;l++){ A=2*l*T+Z; B=A+T; F1=Re(A,amplitude); F2=Im(A,amplitude); X1=FMul(Co ,Re(B,amplitude));
X2=FMul(Co[LS],Im(B,amplitude)); X3=FMul(Co[LS],Re(B,amplitude)); X4=FMul(Co ,Im(B,amplitude));
setRe(A,FSub(FAdd(F1,X1),X2),amplitude); setIm(A,FAdd(FAdd(F2,X3),X4),amplitude); setRe(B,FSub(FAdd(F1,X2),X1),amplitude); setIm(B,FSub(FSub(F2,X3),X4),amplitude); } } T*=2; } } |
||||
| matherp Guru Joined: 11/12/2012 Location: United KingdomPosts: 8584 |
I should probably have mentioned how to do an inverse FFT, i.e. how to convert back from the frequency domain to the amplitude domain. This is really easy using the same CFunction as follows: 1. Swap the real and imaginary components in the frequency domain 2. re-run the FFT 3. Swap the real and imaginary components again 4. Divide all components by the number of elements (2^n) so in our example the code to do this is as follows: OPTION EXPLICIT OPTION DEFAULT NONE dim i% DIM t! DIM amplitudes!(1,7)=(1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0) 'Test data set from rosettacode.org ' ' run the FFT ' FFT 8,amplitudes!() for i%=0 to 7 print "Real=",amplitudes!(0,i%)," Imaginary=",amplitudes!(1,i%) next i% ' ' Now compute the inverse FFT ' for i%=0 to 7 'swap the real and imaginary components t!=amplitudes!(1,i%):amplitudes!(1,i%)=amplitudes!(0,i%):amplitudes!(0,i%)=t! next i% FFT 8,amplitudes!() 'rerun the FFT for i%=0 to 7 'swap the compionents again and divide them by the number of samples t!=amplitudes!(1,i%):amplitudes!(1,i%)=amplitudes!(0,i%):amplitudes!(0,i%)=t! amplitudes!(1,i%)=amplitudes!(1,i%)/8:amplitudes!(0,i%)=amplitudes!(0,i%)/8 print "Real=",amplitudes!(0,i%)," Imaginary=",amplitudes!(1,i%) next i% end The alternative approach to the inverse FFT (just as easy) is to: 1. replace each value in the frequency domain with its complex conjugate (negate the imaginary component) 2. Run the FFT 3. replace each value with its complex conjugate again 4. Divide all components by the number of elements (2^n) This works just as well OPTION EXPLICIT OPTION DEFAULT NONE dim i% DIM t! DIM amplitudes!(1,7)=(1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,0) 'Test data set from rosettacode.org ' ' run the FFT ' FFT 8,amplitudes!() for i%=0 to 7 print "Real=",amplitudes!(0,i%)," Imaginary=",amplitudes!(1,i%) next i% ' ' Now compute the inverse FFT ' for i%=0 to 7 'swap the real and imaginary components amplitudes!(1,i%)= -amplitudes!(1,i%) next i% FFT 8,amplitudes!() 'rerun the FFT for i%=0 to 7 'swap the compionents again and divide them by the number of samples amplitudes!(1,i%)= -amplitudes!(1,i%) amplitudes!(0,i%)=amplitudes!(0,i%)/8 amplitudes!(1,i%)=amplitudes!(1,i%)/8 print "Real=",amplitudes!(0,i%)," Imaginary=",amplitudes!(1,i%) next i% end ' |
||||
| twofingers Guru Joined: 02/06/2014 Location: GermanyPosts: 1133 |
Hi Peter, this is amazing! Thanks for your efforts! 
Michael |
||||
| matherp Guru Joined: 11/12/2012 Location: United KingdomPosts: 8584 |
Oops - left some diagnostic routines in the CFunction - no detrimental impact other than size. Revised version attached: [code] CSUB FFT 'In place FFT for 2^n datasets 00000028 'cos 27BDFFE0 3C039D00 AFB20018 8C62009C 00809021 3C043FC9 AFBF001C AFB10014 AFB00010 24840FDB 8C70006C 0040F809 8C710060 00402021 0220F809 02402821 00402021 0200C821 8FBF001C 8FB20018 8FB10014 8FB00010 03200008 27BD0020 'Re 000420C0 00A42021 03E00008 8C820000 'Im 000420C0 00A42021 03E00008 8C820004 'setRe 000420C0 00C42021 03E00008 AC850000 'setIm 000420C0 00C42021 03E00008 AC850004 'main 27BDFF88 AFB40060 AFB3005C 8C940000 3C139D00 8E620040 AFBF0074 001420C0 AFBE0070 AFB7006C AFB60068 AFB50064 AFB20058 AFB10054 AFB00050 0040F809 00A08021 00408821 00141880 02231821 8E62009C 3C044049 AFA30034 0040F809 24840FDB 00409021 8E620080 00142FC3 02802021 8E770070 8E760064 8E74007C 0040F809 8E75005C 02E0F809 00402021 0040B821 8E62009C 3C044000 0040F809 8E7E0070 03C0F809 00402021 00402821 02C0F809 02E02021 0040B021 3C043EFF 8E62009C 0040F809 3484F2E5 00402821 02A0F809 02C02021 0280F809 00402021 2444FFFF 1880017D AFA20040 00001021 24130002 24420001 0044182A 1460FFFD 00139840 8FA20040 00131843 000227C2 00822021 00042043 00446823 00133083 25A5FFFF AFA30038 18A00170 AFA60018 24020002 AFA20010 00001021 8FA60010 24420001 00063040 0045182A 1460FFFB AFA60010 2484FFFF 18800162 00001021 24160002 24420001 0044182A 1460FFFD 0016B040 19A0001F AE200000 8FA50010 24090001 240C0001 00055FC2 01655821 8FA30010 000B5843 006B5023 014B102A 1440000E 000B1080 01652021 00053080 02201821 02221021 8C680000 00852021 00853823 01284021 0147382A AC480000 00661821 10E0FFF8 00461021 258C0001 01AC102A 14400003 00094840 1000FFE6 01602821 3C149D00 8E82009C 3C044000 8E950064 0040F809 8E970058 00402821 02E0F809 02402021 00409021 8E820080 02602021 0040F809 00132FC3 00402821 02A0F809 02402021 8FA40018 0480002B 0040A821 00009021 3C139D00 8E620080 02402021 00122FC3 0040F809 8E740058 00402021 0280F809 02A02821 00402021 0411FF3B 00000000 8FA50034 00121880 0040A021 00A31821 8E620060 00002021 02802821 0040F809 AC740000 8FA60038 8FA30034 00D22023 00042080 00642021 8E630060 02802821 AC820000 0060F809 00002021 8FA50038 8FA60034 02452021 8FA50018 00042080 26520001 00C42021 00B2182A 1060FFD9 AC820000 0220F021 0000B821 AFA00014 0220A821 10000004 00009021 0256102A 10400034 26B50004 8EB40000 8FC20000 72D41802 02F29821 0062A021 0274102A 1040FFF6 26520001 02802021 02002821 0411FF21 00000000 02602021 02002821 AFA20048 0411FF1C 00000000 00402821 02802021 02003021 0411FF1F 00000000 8FA30048 02003021 00602821 02602021 0411FF19 00000000 02802021 02002821 0411FF11 00000000 02602021 02002821 AFA20048 0411FF0C 00000000 00402821 02802021 02003021 0411FF0F 00000000 8FA30048 02602021 00602821 02003021 0411FF09 00000000 0256102A 1440FFCE 26B50004 8FA40014 8FA50010 24840001 0085102A AFA40014 02F6B821 1440FFC0 27DE0004 8FA60040 18C000AA 8FA30018 8FA50034 00031080 24040001 00A21021 AFA0003C AFA40020 AFA20044 3C119D00 8FA30038 8FA40020 8FA60034 0064001A 008001F4 8FA20044 00042840 AFA50028 AFA60018 AFA2001C AFA0002C 00001812 AFA30024 00031880 AFA30030 8FA40024 1880007B 8FA3002C 8FB2002C 0000A021 02402021 02002821 0411FECF 00000000 8FA60020 02402021 02002821 02469821 AFA20010 0411FECC 00000000 0040B821 8FA20018 02602021 8C560000 02002821 8E350058 0411FEC0 00000000 00402821 02A0F809 02C02021 8FA3001C 02602021 8C760000 02002821 0040F021 8E350058 0411FEB9 00000000 00402821 02A0F809 02C02021 8FA4001C 02002821 8C960000 02602021 8E350058 AFA20014 0411FEAA 00000000 00402821 02A0F809 02C02021 8FA50018 02602021 8CA30000 02002821 8E350058 0040B021 AFA30048 0411FEA2 00000000 8FA30048 00402821 02A0F809 00602021 8E230060 0040A821 8FA40010 8E22005C 03C02821 0040F809 AFA30048 8FA30048 8FA50014 0060F809 00402021 00402821 02003021 02402021 0411FE92 00000000 8E23005C 02E02021 02C02821 0060F809 AFA30048 8FA30048 00402021 0060F809 02A02821 00402821 02003021 02402021 0411FE88 00000000 8E230060 8E22005C 8FA40010 8FA50014 0040F809 AFA30048 8FA30048 03C02821 0060F809 00402021 00402821 02003021 02602021 0411FE75 00000000 8E3E0060 02E02021 03C0F809 02C02821 00402021 03C0F809 02A02821 00402821 02003021 02602021 0411FE6D 00000000 8FA60028 8FA20024 26940001 1682FF8A 02469021 8FA3002C 8FA40020 24630001 8FA50018 8FA60030 AFA3002C 0064102A 8FA3001C 00A62821 00661821 AFA50018 1440FF78 AFA3001C 8FA4003C 8FA50040 24840001 0085102A 10400004 AFA4003C 8FA60028 1000FF60 AFA60020 8FBF0074 8FBE0070 8FB7006C 8FB60068 8FB50064 8FB40060 8FB3005C 8FB20058 8FB10054 8FB00050 03E00008 27BD0078 1000FE89 24130002 1000FEA3 24160002 24020002 1000FE98 AFA20010 End CSUB [/code] #define Version 100 //Version 1.00
#include "../cfunctions.h" float cos(float x){ return FSin(FSub(LoadFloat(0x3fc90fdb),x)); } float Re(int N, float amplitude[]){ return amplitude[N*2]; } float Im(int N, float amplitude[]){ return amplitude[N*2+1]; } void setRe(int N, float val, float amplitude[]){ amplitude[N*2]=val; } void setIm(int N, float val, float amplitude[]){ amplitude[N*2+1]=val; } void main(long *n, float *amplitude){ int P=*n,R1,R2,R3,R4,R,l,j,S1,S2,P1,P2,DV,DP,HA,PT,IP,H,G,T,stage,D,Z,L,LS,A,B; int *V=GetTempMemory(P * 8); float *Co; Co=(void*)V+P*4;; float pi=LoadFloat(0x40490FDB),K,COX,temp,F1,F2,X1,X2,X3,X4; int S=FloatToInt(FAdd(FDiv(FLog(IntToFloat(P)),FLog(LoadFloat(0x40000000))),LoadFloat(0x3EFFF2E5))); R1=2; for(l=0;l<S-1;l++)R1*=2; R=R1-1; R2=R1/2; R4=R1/4; R3=R4+R2; S2=S/2; S1=S-S2; P1=2; for(l=0;l<S1-1;l++)P1*=2; P2=2; for(l=0;l<S2-1;l++)P2*=2; V[0]=0; DV=1; DP=P1; for(j=1;j<=S1;j++){ HA=DP/2; PT=P1-HA; for(l=HA;l<=PT;l+=DP)V =V[l-HA]+DV;
DV*=2; DP=HA; } K=FDiv(FMul(pi,LoadFloat(0x40000000)),IntToFloat(R1)); for(l=0;l<=R4;l++){ COX=cos(FMul(IntToFloat(l),K)); Co =COX;
Co[R2-l]=FSub(0,COX); Co[R2+l]=FSub(0,COX); } for(l=0;l<P1;l++){ IP=l*P2; for(j=0;j<P2;j++){ H=IP+j; G=V[j]*P2+V ;
if(G>H){ temp=Re(G,amplitude); setRe(G,Re(H,amplitude),amplitude); setRe(H,temp,amplitude); temp=Im(G,amplitude); setIm(G,Im(H,amplitude),amplitude); setIm(H,temp,amplitude); } } } T=1; for(stage=0;stage<S;stage++){ D=R2/T; for(Z=0;Z<T;Z++){ L=D*Z; LS=L+R4; for(l=0;l<D;l++){ A=2*l*T+Z; B=A+T; F1=Re(A,amplitude); F2=Im(A,amplitude); X1=FMul(Co ,Re(B,amplitude));
X2=FMul(Co[LS],Im(B,amplitude)); X3=FMul(Co[LS],Re(B,amplitude)); X4=FMul(Co ,Im(B,amplitude));
setRe(A,FSub(FAdd(F1,X1),X2),amplitude); setIm(A,FAdd(FAdd(F2,X3),X4),amplitude); setRe(B,FSub(FAdd(F1,X2),X1),amplitude); setIm(B,FSub(FSub(F2,X3),X4),amplitude); } } T*=2; } } |
||||
| isochronic Guru Joined: 21/01/2012 Location: AustraliaPosts: 689 |
I guess you considered the pic32mx / dspic33 signal-processing library that is available already..What times do you get with say a 512 size dataset ? |
||||
| matherp Guru Joined: 11/12/2012 Location: United KingdomPosts: 8584 |
You can't use libraries in CFunctions MM+, 100MHz, N=512, 29msec |
||||
| isochronic Guru Joined: 21/01/2012 Location: AustraliaPosts: 689 |
Wow...a bit like putting a Concorde engine on a go-cart... [ although I would not use interpreted Basic to manage signal processing in any circumstances !! Why would anyone... just use the C... rant,rant,rant  etc ]
|
||||
 Print this page |
