Eu implementei um radix-4 FFT de 4 pontos e descobri que precisava fazer alguma manipulação dos termos de saída para que ele correspondesse a um dft.
Meu código é uma implementação bastante direta da formulação matricial, então não sei ao certo qual é o problema
// |
// radix-4 butterfly matrix form | complex multiplication
// |
// +- -+ +- -+ | a+ib
// X[0] = | 1 1 1 1 | |x[0]| | * c+id
// X[1] = | 1 -i -1 i | |x[1]| | -------
// X[2] = | 1 -1 1 -1 | |x[2]| | ac + ibc
// X[3] = | 1 i -1 -i | |x[3]| | iad - bd
// +- -+ +- -+ | ------------------
// | (ac-bd) + i(bc+ad)
// |
Alguém pode identificar onde eu errei?
Obrigado,
-David
typedef double fp; // base floating-point type
// naiive N-point DFT implementation as reference to check fft implementation against
//
void dft(int inv, struct cfp *x, struct cfp *y, int N) {
long int i, j;
struct cfp w;
fp ang;
for(i=0; i<N; i++) { // do N-point FFT/IFFT
y[i].r = y[i].i = 0;
if (inv) ang = 2*PI*(fp)i/(fp)N;
else ang = -2*PI*(fp)i/(fp)N;
for (j=0; j<N; j++) {
w.r = cos(j*ang);
w.i = sin(j*ang);
y[i].r += (x[j].r * w.r - x[j].i * w.i);
y[i].i += (x[j].r * w.i + x[j].i * w.r);
}
}
// scale output in the case of an IFFT
if (inv) {
for (i=0; i<N; i++) {
y[i].r = y[i].r/(fp)N;
y[i].i = y[i].i/(fp)N;
}
}
} // dft()
void r4fft4(int inv, int reorder, struct cfp *x, struct cfp *y) {
struct cfp x1[4], w[4];
fp ang, temp;
int i;
// |
// radix-4 butterfly matrix form | complex multiplication
// |
// +- -+ +- -+ | a+ib
// y[0] = | 1 1 1 1 | |x[0]| | * c+id
// y[1] = | 1 -i -1 i | |x[1]| | -------
// y[2] = | 1 -1 1 -1 | |x[2]| | ac + ibc
// y[3] = | 1 i -1 -i | |x[3]| | iad - bd
// +- -+ +- -+ | ------------------
// | (ac-bd) + i(bc+ad)
// |
if (inv) ang = 2*PI/(fp)4; // invert sign for IFFT
else ang = -2*PI/(fp)4;
//
w[1].r = cos(ang*1); w[1].i = sin(ang*1); // twiddle1 = exp(-2*pi/4 * 1);
w[2].r = cos(ang*2); w[2].i = sin(ang*2); // twiddle2 = exp(-2*pi/4 * 2);
w[3].r = cos(ang*3); w[3].i = sin(ang*3); // twiddle3 = exp(-2*pi/4 * 3);
// *1 *1 *1 *1
y[0].r = x[0].r + x[1].r + x[2].r + x[3].r;
y[0].i = x[0].i + x[1].i + x[2].i + x[3].i;
// *1 *-i *-1 *i
x1[1].r = x[0].r + x[1].i - x[2].r - x[3].i;
x1[1].i = x[0].i - x[1].r - x[2].i + x[3].r;
// *1 *-1 *1 *-1
x1[2].r = x[0].r - x[1].r + x[2].r - x[3].r;
x1[2].i = x[0].i - x[1].i + x[2].i - x[3].i;
// *1 *i *-1 *-i
x1[3].r = x[0].r - x[1].i - x[2].r + x[3].i;
x1[3].i = x[0].i + x[1].r - x[2].i - x[3].r;
//
y[1].r = x1[1].r*w[1].r - x1[1].i*w[1].i; // scale radix-4 output
y[1].i = x1[1].i*w[1].r + x1[1].r*w[1].i;
//
y[2].r = x1[2].r*w[2].r - x1[2].i*w[2].i; // scale radix-4 output
y[2].i = x1[2].i*w[2].r + x1[2].r*w[2].i;
//
y[3].r = x1[3].r*w[3].r - x1[3].i*w[3].i; // scale radix-4 output
y[3].i = x1[3].i*w[3].r + x1[3].r*w[3].i;
// reorder output stage ... mystery as to why I need this
if (reorder) {
temp = y[1].r;
y[1].r = -1*y[1].i;
y[1].i = temp;
//
y[2].r = -1*y[2].r;
//
temp = y[3].r;
y[3].r = y[3].i;
y[3].i = -1*temp;
}
// scale output for inverse FFT
if (inv) {
for (i=0; i<4; i++) { // scale output by 1/N for IFFT
y[i].r = y[i].r/(fp)4;
y[i].i = y[i].i/(fp)4;
}
}
} // r4fft4()
Respostas:
Acabei de portar um radix-4 DIF fft do código S. Burrus Fortran para Java. Na verdade, carece de várias otimizações, antes de tudo o fator de twiddle acionado por tabela (os fatores sin e cos devem ser pré-calculados). Isso deve acelerar o fft um pouco mais (talvez 50%). Eu tenho que cortar um pouco para isso, mas se alguém tiver a resposta correta, ficarei muito feliz e agradecido. Vou postar o código otimizado o mais cedo possível Espero que talvez com alguns testes de velocidade vs algoritmo radix-2.
Além disso, as multiplicações por 1 e sqrt (-1) não são removidas. Removê-los irá acelerar um pouco mais. Mas, em geral, o IMHO radix-4 parece não ser mais de 25% mais rápido que um radix-2, então não sei se a relação velocidade / complexidade vale realmente a pena. Lembre-se de que bibliotecas muito otimizadas, como a FFTW, estão amplamente disponíveis e usadas, portanto esse esforço pode ser apenas um "desvio" pessoal!
Aqui está o código java. Portá-lo para C, C ++ ou C # deve ser muito fácil.
Aqui está o código Radix-4 DIF FORTRAN original de Sidney Burrus:
Radix-4, DIF, uma borboleta FFT
fonte
Primeiro, sua suposta 'borboleta radix-4' é uma DFT de 4 pontos, não uma FFT. Possui 16 operações complexas (ie: N ao quadrado). Uma FFT típica de 4 pontos teria apenas Nlog (base 2) N (= 8 para N = 4). Segundo, você tem alguns supostos fatores w [] .r e w [] .i 'scale' que não pertencem. Talvez você os tenha obtido a partir de uma borboleta radix-4 mostrada em um gráfico maior. Uma borboleta assim teria alguns acréscimos nos estágios anexados, mas na verdade não fazem parte da borboleta. Uma FFT de 4 pontos possui apenas uma borboleta interna de -j quando projetada para uma FFT de expoente negativo.
Em vez de tentar consertar seu código, é tão fácil escrever o meu, como mostrado abaixo (compilador DevC ++; saídas anexadas no final do código):
Primeiro, os dados de entrada (4 reais, 4 imaginários) são impressos. Em seguida, é tomada uma DFT de 4 pontos. Os resultados (ano [] e ano [] mais amp / fase) são impressos. Como os dados originais r [] e i [] não foram sobrescritos ao fazer a DFT, essas entradas são reutilizadas como entradas para a FFT de 4 pontos. Observe que o último possui menos operações +/- do que o DFT.
O código para a FFT não é particularmente elegante nem eficiente - existem muitas maneiras de fazer borboletas. O código acima corresponde às quatro borboletas radix-2 mostradas no livro de Rabiner e Gold “Teoria e Aplicação do Processamento Digital de Sinais” (p. 580, Fig. 10.9), com twiddles modificados para refletir um expoente negativo (os usados para o número do livro foram positivos). Observe que há apenas um toque de -j no código, e isso não requer uma multiplicação (é uma troca de troca / sinal).
Após a FFT, os resultados são impressos. Eles são os mesmos que o DFT
E, finalmente, os resultados em escala da FFT são usados como entradas para uma FFT inversa. Isso é realizado através do método de 'troca' ou 'reversão da lista' (ou seja: se FFT (r, i) é uma FFT direta, então a FFT (i, r) é inversa - desde que, naturalmente, a FFT seja capaz de lidar com entradas / saídas complexas - em outras palavras - sem rotinas 'apenas reais', que geralmente presumem que as entradas imaginárias são zero). Este método foi descrito há quase 25 anos em:
P. Duhamel, B. Piron, JM Etcheto, “On Computing the DFT Inverso”, IEEE Transactions on Acoustics, Speech and Signal Processing, vol. 36, fevereiro de 1988, pp. 285-286.
O resultado do inverso é então impresso. É o mesmo que os dados de entrada originais.
fonte