DSP - Filtrado en el dominio de la frecuencia a través de FFT

Question

He estado jugando un poco con la implementación Exocortex de la FFT, pero estoy teniendo algunos problemas.

Cada vez que modifico las amplitudes de los bins de frecuencia antes de llamar al iFFT, la señal resultante contiene algunos clics y pops, especialmente cuando hay bajas frecuencias en la señal (como baterías o bajos). Sin embargo, esto no sucede si atenúo todos los contenedores por el mismo factor.

me deja poner un ejemplo de la memoria intermedia de salida de un FFT 4-muestra:

// Bin 0 (DC) 
FFTOut[0] = 0.0000610351563 
FFTOut[1] = 0.0 

// Bin 1 
FFTOut[2] = 0.000331878662 
FFTOut[3] = 0.000629425049 

// Bin 2 
FFTOut[4] = -0.0000381469727 
FFTOut[5] = 0.0 

// Bin 3, this is the first and only negative frequency bin. 
FFTOut[6] = 0.000331878662 
FFTOut[7] = -0.000629425049 

La salida se compone de pares de flotadores, cada uno representando las partes real e imaginay de un solo bin. Entonces, bin 0 (índices de matriz 0, 1) representaría las partes real e imaginaria de la frecuencia de DC. Como puede ver, los contenedores 1 y 3 tienen los mismos valores (excepto el signo de la parte Im), así que supongo que bin 3 es la primera frecuencia negativa, y finalmente los índices (4, 5) serían los últimos positivos compartimiento de frecuencia.

A continuación, para atenuar el contenedor de frecuencia 1 esto es lo que hago:

// Attenuate the 'positive' bin 
FFTOut[2] *= 0.5; 
FFTOut[3] *= 0.5; 

// Attenuate its corresponding negative bin. 
FFTOut[6] *= 0.5; 
FFTOut[7] *= 0.5; 

Para las pruebas reales que estoy usando una FFT de 1024 de longitud y siempre proporcionan todas las muestras por lo que no es 0-padding necesario.

// Attenuate 
var halfSize = fftWindowLength/2; 
float leftFreq = 0f; 
float rightFreq = 22050f; 
for(var c = 1; c < halfSize; c++) 
{ 
    var freq = c * (44100d/halfSize); 

    // Calc. positive and negative frequency indexes. 
    var k = c * 2; 
    var nk = (fftWindowLength - c) * 2; 

    // This kind of attenuation corresponds to a high-pass filter. 
    // The attenuation at the transition band is linearly applied, could 
    // this be the cause of the distortion of low frequencies? 
    var attn = (freq < leftFreq) ? 
        0 : 
        (freq < rightFreq) ? 
         ((freq - leftFreq)/(rightFreq - leftFreq)) : 
         1; 

    // Attenuate positive and negative bins. 
    mFFTOut[ k ] *= (float)attn; 
    mFFTOut[ k + 1 ] *= (float)attn; 
    mFFTOut[ nk ] *= (float)attn; 
    mFFTOut[ nk + 1 ] *= (float)attn; 
} 

Obviamente estoy haciendo algo mal pero no puedo entender qué.

No quiero usar la salida FFT como un medio para generar un conjunto de coeficientes FIR ya que estoy tratando de implementar un ecualizador dinámico muy básico.

¿Cuál es la forma correcta de filtrar en el dominio de la frecuencia? ¿Qué me estoy perdiendo?

Además, ¿es realmente necesario atenuar las frecuencias negativas también? He visto una implementación de FFT donde neg. los valores de frecuencia se ponen a cero antes de la síntesis.

Gracias de antemano.

Answer 1

Hay dos problemas: la forma de usar la FFT y el filtro en particular.

El filtrado se implementa tradicionalmente como convolución en el dominio del tiempo. Tienes razón en que multiplicar los espectros de las señales de entrada y filtro es equivalente. Sin embargo, cuando utiliza la Transformada de Fourier Discreta (DFT) (implementada con un algoritmo de Transformada Rápida de Fourier para velocidad), realmente calcula una versión muestreada del espectro verdadero. Esto tiene muchas implicaciones, pero la más relevante para el filtrado es la implicación de que la señal del dominio del tiempo es periódica.

Aquí hay un ejemplo. Considere una señal de entrada sinusoidal x con 1.5 ciclos en el período, y un filtro de paso bajo simple h. En Matlab sintaxis/Octava:

N = 1024; 
n = (1:N)'-1; %'# define the time index 
x = sin(2*pi*1.5*n/N); %# input with 1.5 cycles per 1024 points 
h = hanning(129) .* sinc(0.25*(-64:1:64)'); %'# windowed sinc LPF, Fc = pi/4 
h = [h./sum(h)]; %# normalize DC gain 

y = ifft(fft(x) .* fft(h,N)); %# inverse FT of product of sampled spectra 
y = real(y); %# due to numerical error, y has a tiny imaginary part 
%# Depending on your FT/IFT implementation, might have to scale by N or 1/N here 
plot(y); 

Y aquí está el gráfico:

El problema técnico al inicio del bloque no es lo que esperamos en absoluto. Pero si considera fft(x), tiene sentido. La Transformada de Fourier Discreta supone que la señal es periódica dentro del bloque de transformación. Por lo que sabe la DFT, nos preguntamos por la transformación de un periodo de esto:

Esto lleva a la primera consideración importante cuando se filtra con DFT: en realidad se está implementando circular convolution, no convolución lineal. Entonces, el "error" en el primer gráfico no es realmente un problema cuando se consideran las matemáticas. Entonces, la pregunta es: ¿hay alguna manera de evitar la periodicidad? La respuesta es sí: use overlap-save processing. Básicamente, calcula los productos N-long como arriba, pero solo mantiene N/2 puntos.

Nproc = 512; 
xproc = zeros(2*Nproc,1); %# initialize temp buffer 
idx = 1:Nproc; %# initialize half-buffer index 
ycorrect = zeros(2*Nproc,1); %# initialize destination 
for ctr = 1:(length(x)/Nproc) %# iterate over x 512 points at a time 
    xproc(1:Nproc) = xproc((Nproc+1):end); %# shift 2nd half of last iteration to 1st half of this iteration 
    xproc((Nproc+1):end) = x(idx); %# fill 2nd half of this iteration with new data 
    yproc = ifft(fft(xproc) .* fft(h,2*Nproc)); %# calculate new buffer 
    ycorrect(idx) = real(yproc((Nproc+1):end)); %# keep 2nd half of new buffer 
    idx = idx + Nproc; %# step half-buffer index 
end 

Y aquí está la gráfica de ycorrect:

Este cuadro tiene sentido - esperamos un transitorio de arranque desde el filtro, entonces el resultado se asienta en la respuesta sinusoidal de estado estacionario. Tenga en cuenta que ahora x puede ser arbitrariamente largo. La limitación es Nproc > 2*min(length(x),length(h)).

Ahora en el segundo problema: el filtro en particular.En el bucle, se crea un filtro que es el espectro es esencialmente H = [0 (1:511)/512 1 (511:-1:1)/512]'; Si lo hace hraw = real(ifft(H)); plot(hraw), se obtiene:

Es difícil de ver, pero hay un montón de no-cero puntos en el borde extremo izquierdo de la gráfica , y luego un grupo más en el extremo derecho. Utilizando una función de freqz función de Octave para mirar la respuesta de frecuencia vemos (haciendo freqz(hraw)):

La respuesta de magnitud tiene una gran cantidad de ondas de la envoltura de paso alto a cero. Nuevamente, la periodicidad inherente en el DFT está en acción. En lo que respecta al DFT, hraw se repite una y otra vez. Pero si toma un período de hraw, como freqz, su espectro es bastante diferente de la versión periódica.

Definamos una nueva señal: hrot = [hraw(513:end) ; hraw(1:512)]; Simplemente giramos la salida de DFT sin procesar para que sea continua dentro del bloque. Ahora veamos la respuesta de frecuencia usando freqz(hrot):

Mucho mejor. El sobre deseado está allí, sin todas las ondas. Por supuesto, la implementación no es tan simple ahora, tiene que hacer una multiplicación compleja completa por fft(hrot) en lugar de simplemente escalar cada contenedor complejo, pero al menos obtendrá la respuesta correcta.

Tenga en cuenta que para la velocidad, normalmente precalcular el DFT del h acolchado, lo dejé solo en el circuito para compararlo más fácilmente con el original.

Answer 2

Primero, sobre la normalización: es un problema conocido (no). El DFT/IDFT requeriría un factor 1/sqrt (N) (aparte de los factores coseno/seno estándar) en cada uno (directo y inverso) para hacerlos simétricos y verdaderamente invertibles. Otra posibilidad es dividir uno de ellos (el directo o el inverso) por N, esto es una cuestión de conveniencia y gusto. A menudo las rutinas de FFT no realizan esta normalización, se espera que el usuario lo conozca y se normalice como lo prefiera. See

Segundo: en un (digamos) de 16 puntos DFT, lo que se llama a la bin 0 correspondería a la frecuencia cero (DC), bin 1 baja frec ... bin 4 frec medio, bin 8 a la frecuencia más alta y bins 9 ... 15 a las "frecuencias negativas". En su ejemplo, entonces, bin 1 es en realidad la frecuencia baja y la frecuencia media. Aparte de esta consideración, no hay nada conceptualmente incorrecto en su "igualación". No entiendo lo que quiere decir con "la señal se distorsiona a bajas frecuencias". ¿Cómo observa eso?

Answer 3

¿Está atenuando el valor de la muestra de frecuencia de CC a cero? Parece que no estás atenuándolo en absoluto en tu ejemplo. Como está implementando un filtro de paso alto, también debe establecer el valor de CC en cero.

Esto explicaría la distorsión de baja frecuencia. Tendría muchas ondulaciones en la respuesta de frecuencia a bajas frecuencias si ese valor de CC no es cero debido a la gran transición.

Aquí se muestra un ejemplo en MATLAB/Octave para demostrar lo que podría estar sucediendo:

N = 32; 
os = 4; 
Fs = 1000; 
X = [ones(1,4) linspace(1,0,8) zeros(1,3) 1 zeros(1,4) linspace(0,1,8) ones(1,4)]; 
x = ifftshift(ifft(X)); 
Xos = fft(x, N*os); 
f1 = linspace(-Fs/2, Fs/2-Fs/N, N); 
f2 = linspace(-Fs/2, Fs/2-Fs/(N*os), N*os); 

hold off; 
plot(f2, abs(Xos), '-o'); 
hold on; 
grid on; 
plot(f1, abs(X), '-ro'); 
hold off; 
xlabel('Frequency (Hz)'); 
ylabel('Magnitude'); 

frequency response http://www.freeimagehosting.net/uploads/e10109e535.png

en cuenta que en mi código, estoy creando un ejemplo del valor de corriente continua siendo distinto de cero , seguido por un cambio abrupto a cero, y luego una aceleración. Luego tomo el IFFT para transformarme en el dominio del tiempo. Luego realizo un fft de relleno cero (que MATLAB realiza automáticamente cuando pasas en un tamaño de fft mayor que la señal de entrada) en esa señal de dominio de tiempo. El relleno cero en el dominio del tiempo da como resultado la interpolación en el dominio de la frecuencia. Al usar esto, podemos ver cómo responderá el filtro entre muestras de filtros.

Una de las cosas más importantes para recordar es que, aunque está configurando valores de respuesta de filtro en frecuencias determinadas al atenuar las salidas de la DFT, esto no garantiza nada para las frecuencias que ocurren entre puntos de muestra. Esto significa que cuanto más abruptos sean sus cambios, mayor será el sobreimpulso y la oscilación entre las muestras.

Ahora para responder a su pregunta sobre cómo se debe hacer este filtrado. Hay varias formas, pero una de las más fáciles de implementar y comprender es el método de diseño de ventanas. El problema con su diseño actual es que el ancho de transición es enorme. La mayoría de las veces, deseará las transiciones más rápidas posible, con la menor onda posible.

En el siguiente código, voy a crear un filtro ideal y mostrar la respuesta:

N = 32; 
os = 4; 
Fs = 1000; 
X = [ones(1,8) zeros(1,16) ones(1,8)]; 
x = ifftshift(ifft(X)); 
Xos = fft(x, N*os); 
f1 = linspace(-Fs/2, Fs/2-Fs/N, N); 
f2 = linspace(-Fs/2, Fs/2-Fs/(N*os), N*os); 

hold off; 
plot(f2, abs(Xos), '-o'); 
hold on; 
grid on; 
plot(f1, abs(X), '-ro'); 
hold off; 
xlabel('Frequency (Hz)'); 
ylabel('Magnitude'); 

frequency response http://www.freeimagehosting.net/uploads/c86f5f1700.png

en cuenta que hay una gran cantidad de oscilación causada por los cambios bruscos.

La Transformada de Fourier Discreta o FFT es una versión de muestra de la Transformada de Fourier. La Transformada de Fourier se aplica a una señal en el rango continuo -infinito hasta el infinito, mientras que el DFT se aplica sobre un número finito de muestras. Esto en efecto da como resultado una ventana cuadrada (truncamiento) en el dominio de tiempo cuando se usa DFT ya que solo estamos tratando con un número finito de muestras. Desafortunadamente, el DFT de una onda cuadrada es una función de tipo sinc (sin (x)/x).

El problema de tener transiciones nítidas en el filtro (salto rápido de 0 a 1 en una muestra) es que tiene una respuesta muy larga en el dominio de tiempo, que se trunca con una ventana cuadrada. Para ayudar a minimizar este problema, podemos multiplicar la señal del dominio del tiempo por una ventana más gradual.Si multiplicamos una ventana Hanning añadiendo la línea:

x = x .* hanning(1,N).'; 

después de tomar la IFFT, obtenemos esta respuesta:

frequency response http://www.freeimagehosting.net/uploads/944da9dd93.png

por lo que recomiendo tratando de poner en práctica el método de diseño de la ventana, ya que es bastante simple (hay mejores formas, pero se vuelven más complicadas). Como está implementando un ecualizador, supongo que desea poder cambiar las atenuaciones sobre la marcha, por lo que le sugiero que calcule y almacene el filtro en el dominio de la frecuencia siempre que haya un cambio en los parámetros, y luego puede aplicarlo a cada búfer de audio de entrada tomando el fft del búfer de entrada, multiplicando por sus muestras de filtro de dominio de frecuencia y luego realizando el ifft para regresar al dominio de tiempo. Esto será mucho más eficiente que todas las ramificaciones que está haciendo para cada muestra.

Answer 4

Su problema principal es que las frecuencias no están bien definidas en intervalos de tiempo cortos. Esto es particularmente cierto para las frecuencias bajas, que es la razón por la que se nota el problema más allí.

Por lo tanto, cuando saca segmentos realmente cortos del tren de sonido y los filtra, los segmentos filtrados no se filtran de forma que produzca una forma de onda continua, y escucha los saltos entre segmentos y esto es lo que genera los clics aquí.

Por ejemplo, tomar algunos números razonables: Comienzo con una forma de onda a 27,5 Hz (A0 en un piano), se digitalizaron a 44.100 Hz, que se verá como este (donde la parte roja es de 1024 muestras de longitud):

alt text http://i48.tinypic.com/zim802.png

Así que primero comenzaremos con un paso bajo de 40Hz. Por lo tanto, dado que la frecuencia original es inferior a 40 Hz, un filtro de paso bajo con un corte de 40 Hz no debería tener ningún efecto, y obtendremos una salida que coincide casi exactamente con la entrada. ¿Derecha? Incorrecto, incorrecto, incorrecto - y este es básicamente el núcleo de su problema. El problema es que para las secciones cortas, la idea de 27.5 Hz no está claramente definida, y no se puede representar bien en el DFT.

Eso 27,5 Hz no es particularmente significativo en el segmento corto se puede ver mirando el DFT en la figura a continuación. Tenga en cuenta que aunque el DFT del segmento más largo (puntos negros) muestra un pico a 27.5 Hz, el corto (puntos rojos) no lo hace.

alt text http://i50.tinypic.com/14w6luw.png

Claramente, a continuación, filtrar a continuación 40Hz, se acaba de capturar el desplazamiento de CC, y el resultado del filtro de paso bajo de 40 Hz se muestra en verde a continuación.

alt text http://i48.tinypic.com/2vao21w.png

La curva azul (tomado con un 200 Hz de corte) está empezando a coincidir mucho mejor. Pero tenga en cuenta que no son las bajas frecuencias las que lo hacen coincidir bien, sino la inclusión de altas frecuencias. No es hasta que incluimos todas las frecuencias posibles en el segmento corto, hasta 22 KHz, que finalmente obtenemos una buena representación de la onda sinusoidal original.

La razón de todo esto es que un pequeño segmento de una onda sinusoidal de 27,5 Hz es no una onda senoidal de 27,5 Hz, y su DFT no tiene mucho que ver con 27,5 Hz.