Optimización de algoritmos basados en FFT Si necesita un programa completo en lenguaje C (utilizando las propiedades de los factores de rotación), deje un mensaje, ¡gracias! ! ! (Hay un código central, por favor dame algún consejo)
Implementación (descripción C)
#include lt;stdio.hgt;
#include lt;math.hgt;
# incluir lt;stdlib.hgt;
//#include "complex.h"
// ------------------ -------------------------------------------------- ------
#define N 8 //64
#define M 3 //6 //2^m=N
#define PI 3.1415926
// --------------------------------------- ----------------------------------
giro flotante[N/2] = { 1,0, 0,707, 0,0, -0,707};
flotante x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0}
flotante x_i; [N]; //N=8
/*
giro flotante[N/2] = {1, 0,9951, 0,9808, 0,9570, 0,9239, 0,8820, 0,8317, 0,7733,
0.7075, 0.6349, 0.5561, 0.4721, 0.3835, 0.2912, 0.1961, 0.0991,
0.0000, -0.0991, -0.1961, -0.2912, -0.3835, -0. 1,- 0,5561, -0,6349,
-0,7075, -0,7733, 0,8317, -0,8820, -0,9239, -0,9570, -0,9808, -0,9951} //N=64
flotante x_r [N]={1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1,
0, 0 , 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0 , 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,};
flotante x_i[N];
*/
ARCHIVO *fp ;
// ---------------------------------- función ---- - ---------------------------------
/**
* Inicialización Salida parte imaginaria
*/
static void fft_init( void )
{
int i;
for(i=0; ilt; N; i) x_i[i] = 0.0;
}
/**
* Invertir el algoritmo
* Este algoritmo tiene margen de mejora
*/
bitrev vacío estático (void)
{.
int p=1, q, i;
int bit_rev[ N ]; //
float xx_r[ N ];
bit_rev[ 0 ] = 0;
mientras( p lt; N )
{
for(q=0; qlt; p; q )
{
bit_rev[ q ] = bit_rev[ q ] * 2;
bit_rev[ q p ] = bit_rev[ q ] 1;
}
p *= 2;
}
for(i=0; ilt; N; i) xx_r[ i ] = x_r[ i ];
for(i=0; ilt; N; i ) x_r[i] = xx_r[ bit_rev[i] ]; /* ------------ agregar por sshc625 ------------ */
static void bitrev2( void )
{
return ;
}
/* */
visualización vacía( void )
{
printf("\n\n");
int i
for(i=0; ilt; N; i)
printf("f\tf\n", x_r[i], x_i[i]);
}
/**
*
*/
void fft1( void )
{ fp = fopen("log1.txt", "a "); > int L, i, b, j, p, k, tx1, tx2;
float TR, TI, temp;
// Variables temporales
float tw1, tw2;
/* Deep M. Recorre las capas L es la capa actual y el número total de capas es M. */
for(L=1; Llt;=M; L)
{
fprintf(fp, "---------- Layer=d- ---------\n", L);
/* b es de gran importancia, b representa el número de puntos de muestra de entrada que las partículas de la capa actual tener*/
b = 1;
i = L - 1;
mientras(i gt; 0)
{
b *= 2;
i--;
}
// ----------- --- Si la capa exterior hace circular las partículas, la capa interior es más lógicamente circular para los puntos de muestra --------------
/*
* outter es más lógico para participar en DFT Los puntos de muestra se ciclan
* L=1, se cicla 1 vez (4 partículas, 2 puntos de muestra por partícula)
* L= 2, ciclo 2 veces (2 partículas, 4 puntos de muestra por partícula)
* L=3, ciclo 4 veces (1 partícula, 8 puntos de muestra por partícula)
*/
for(j=0; jlt; b; j)
{
/* Encuentra el factor de rotación tw1 */
p = 1;
i = M - L; // M es el número total de capas, L es la capa actual.
while(i gt; 0)
{
p = p*2
i--;
}
p = p * j; /p >
tx1 = p N;
tx2 = tx1 3*N/4;
tx2 = tx2 N
// tw1 es la parte cos, parte real; tw2 es la parte sin, la parte imaginaria.
tw1 = (tx1gt;=N/2)? -twiddle[tx1-N/2]: twiddle[tx1];
tw2 = (tx2gt;=N/2)? -twiddle[tx2-(N/2)]: twiddle[tx2];
/*
* par interno de partículas Loop
* L=1, bucle 4 veces (4 partículas, 2 entradas por partícula)
* L=2, bucle 2 veces (2 partículas, 4 entradas por partícula)
* L=3, bucle una vez (1 partícula, 8 entradas por partícula)
*/
for(k=j; klt; N; k=k 2*b)
{
TR = x_r[k] // TR es A, x_r[k b] es B.
TI = x_i[ k]; p>
temp = x_r[k b];
/*
* Si repasas los dos números complejos (a j*b) (c j*d) ¿Cuáles son los reales? y partes imaginarias después de la multiplicación
* Puedes entender por qué se realiza la siguiente operación La entrada es solo un número real cuando L = 1, y las capas posteriores
* Las entradas son. todos los números complejos Para que las entradas de todas las capas sean números complejos, tenemos que dejar L = 1 cuando
* La parte imaginaria de entrada es 0
* x_i[k b. ]* tw2 es la multiplicación de dos números imaginarios
*/
fprintf(fp, "tw1=f, tw2=f\n", tw1, tw2);
x_r[k] = TR x_r[k b]*tw1 x_i[k b]*tw2;
x_i[k] = TI - x_r[k b]*tw2 x_i[k b]*tw1;
x_r[k b] = TR - x_r[k b]*tw1 - x_i[k b]*tw2;
x_i[k b] = TI temp*tw2 - x_i[k b]* tw1;
fprintf(fp, "k=d, x_r[k]=f, x_i[k]=f\n", k, x_r[k], x_i[k]);
fprintf(fp, "k=d, x_r[k]=f, x_i[k]=f\n", k b, x_r[k b], x_i[k b]);
} / /
} //
} //
}
/**
* --- --------- agregar por sshc625 ------------
* El proceso de implementación es
* for(Layer )
* para (gránulo)
* para (muestra)
*
*
*
*
*/
void fft2( void )
{ fp = fopen("log2.txt", "a ") ;
int cur_layer, gr_num, i, k, p
float tmp_real, tmp_imag, temp; // Variable temporal, registra la parte real
float tw1, tw2; // Factor de rotación, tw1 es la parte real cos parte del factor de rotación, tw2 es la parte imaginaria pecado del factor de rotación Part.
int step; // Paso
int sample_num // El número total de muestras de partículas (diferentes en cada capa, debido a la entrada de partículas) en cada capa es diferente)
/* Recorrer capas*/
for(cur_layer=1; cur_layerlt;=M; cur_layer)
{
/* Encuentra cuántas partículas de la capa actual (gr_num) */
gr_num = 1;
i = M - cur_layer
mientras(i gt; 0)
{
i--;
gr_num *= 2;
/* Muestra de entrada para cada partícula Número N' */
sample_num = (int)pow(2, cur_layer);
/* El paso es N'. /2 */
paso = sample_num/2;
/* */
k = 0
/* Bucle el partículas*/
for(i=0; ilt; gr_num; i )
{
/*
* Realizar un bucle en la muestra puntos, preste atención al límite superior y al paso
*/
for(p=0; plt; sample_num/2; p )
{
// Factor de rotación, necesita optimización...
tw1 = cos(2*PI*p/pow(2, cur_layer));
tw2 = - sin(2*PI*p/pow( 2. cur_layer));
tmp_real = x_r[k p];
tmp_imag = x_i[k p];
temp = x_r[k p paso];
/*(tw1 jtw2)(x_r[k] jx_i[k])
*
* real: tw1*x_r[k] - tw2*x_i [k]
* imag: tw1*x_i[k] tw2*x_r[k]
* No quiero hacer abstracción uno
* typedef struct {
* double real; // parte real
* double imag; // parte imaginaria
* } operaciones complejas y <; /p>
* ¿Se debe a consideraciones de eficiencia simplificar operaciones con números complejos?
*/
/* Algoritmo de mariposa*/
x_r[ k p] = tmp_real ( tw1*x_r[k p paso] - tw2*x_i[k p paso] );
x_i[k p] = tmp_imag ( tw2*x_r[k p paso] tw1*x_i[k p paso] );
/* X[k] = A(k) WB(k)
* X[k N/2] = A(k)-WB(k ) propiedades se puede optimizar aquí*/
// El factor de rotación debe optimizarse...
tw1 = cos(2*PI*(p paso)/pow(2, cur_layer ) );
tw2 = -sin(2*PI*(p paso)/pow(2, cur_layer));
x_r[k p paso] = tmp_real ( tw1*temp - tw2*x_i[k p paso] );
x_i[k p paso] = tmp_imag ( tw2*temp tw1*x_i[k p paso] );
printf("k=d , x_r[k]=f, x_i[k]=f\n", k p, x_r[k p], x_i[k p]);
printf("k=d, x_r[k]= f ， */
k = 2*paso
}
}
}
/*
/*
p>
* Posdata:
* Si las partículas circulan en la capa exterior o las muestras se introducen en la capa exterior, parece similar, la complejidad es exactamente la misma.
* Pero con mis calificaciones Me tomó mucho tiempo descifrar estas docenas de líneas de código.
* Encontré una lección muy útil. Escribí algunos algoritmos hace mucho tiempo, la mayoría de los cuales eran recursivos.
* Reduzca la cantidad de datos, reduzca y reduzca nuevamente, y use métodos inductivos para encontrar. Descubra las reglas para aumentar la cantidad de datos
* Por ejemplo, FFT
* 1. Primero escriba el código de LayerI y luego use la salida de LayerI como entrada de LayerII; y luego escribirlo de nuevo
Código; ......
* Se necesitan alrededor de 3 niveles para calcular estadísticamente el patrón. Esto es lo mismo que la recursividad, ¡simplemente escriba uno o dos niveles primero y saldrá naturalmente!
* 2. Algunas funciones se pueden escribir en pseudocódigo. No se apresure a encontrar el resultado, reduzca la complejidad y agréguelo después de determinar el resultado lógico.
* Para. Por ejemplo, el factor de rotación se puede codificar, simplemente escriba 1.0. Escriba el factor de rotación una vez completado el proceso.
* Solo unas pocas palabras, ¡realmente sudé mucho!
*/
void dft ( void )
{
int i, n, k, tx1, tx2
float tw1, tw2;
float xx_r [N], xx_i[N];
/*
* borra cualquier dato en las matrices de resultados reales e imaginarios antes a DFT
*/
for(k=0; klt;=N-1; k)
xx_r[k] = xx_i[k] = x_i[k] = 0.0;
// calcula el DFT
for(k=0; klt;=(N-1); k )
{
for(n=0; nlt;=(N-1); n )
{
tx1 = (n*k);
tx2 = tx1 (3*N)/ 4;
tx1 = tx1(N
tx2 = tx2(N); > if(tx1 gt; = (N/2))
tw1 = -twiddle[tx1-(N/2)]
else
tw1; = twiddle[tx1];
if (tx2 gt; = (N/2))
tw2 = -twiddle[tx2-(N/2)];
else
tw2 = giro [tx2];
xx_r[k] = xx_r[k] x_r[n]*tw1
xx_i[ k] = xx_i[k] x_r[n]*tw2;
}
xx_i[k] = -xx_i[k]
}
// mostrar
for(i=0; ilt; N; i)
printf("f\tf\n", xx_r[i], xx_i [yo]);
}
// ---------------------------- ----------------- --
----------------------------------
int principal( vacío )
{
fft_init( );
bitrev(
// bitrev2(
//fft1();
fft2( );
mostrar(
sistema( "pausa"
// dft()); ;
p>
return 1;
}
Este artículo proviene del blog de CSDN, indique la fuente al reimprimir: /sshcx/archive/ 2007/06/14/1651616.aspx