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

procesamiento de imagenes con fourier

Question

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

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

Permítame poner un ejemplo del búfer de salida de una FFT de 4 muestras:

// 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 está compuesta por pares de flotadores, cada uno de los cuales representa las partes reales e imaginativas de un solo contenedor. Por lo tanto, bin 0 (índices de matriz 0, 1) representaría las partes reales e imaginarias de la frecuencia de CC. Como puede ver, los intervalos 1 y 3 tienen los mismos valores, (excepto por el signo de la parte Im), así que supongo que el intervalo 3 es la primera frecuencia negativa, y finalmente los índices (4, 5) serían el último positivo bandeja de frecuencia.

Luego para atenuar la frecuencia bin 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, estoy usando una FFT de 1024 de longitud y siempre proporciono todas las muestras, por lo que no se necesita relleno de 0.

// 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 sé 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? lo que me falta

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 por adelantado.

Answer 1

¿Está atenuando el valor de la muestra de frecuencia de CC a cero? Parece que no lo estás atenuando 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 una gran cantidad de rizado en la respuesta de frecuencia a bajas frecuencias si ese valor de CC no fuera cero debido a la gran transición.

Aquí hay 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'');

respuesta de frecuencia http://www.freeimagehosting.net/uploads/e10109e535.png

Observe que en mi código, estoy creando un ejemplo de que el valor de CC no es cero, seguido de un cambio abrupto a cero, y luego una rampa ascendente. Luego tomo el IFFT para transformarme en el dominio del tiempo. Luego realizo un fft sin relleno (que MATLAB realiza automáticamente cuando se pasa en un tamaño de pie más grande 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. Usando esto, podemos ver cómo responderá el filtro entre las muestras de filtro.

Una de las cosas más importantes que debe 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 se producen entre los puntos de muestra. Esto significa que cuanto más abruptos sean sus cambios, más se producirá un exceso y una 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 fluctuación posible.

En el siguiente código, 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''); 

respuesta de frecuencia http://www.freeimagehosting.net/uploads/c86f5f1700.png

Observe que hay mucha oscilación causada por los cambios bruscos.

La FFT o la transformada discreta de Fourier es una versión muestreada de la transformada de Fourier. La transformada de Fourier se aplica a una señal en el rango continuo de infinito a infinito, mientras que la DFT se aplica a un número finito de muestras. En efecto, esto resulta en una ventana cuadrada (truncamiento) en el dominio del tiempo cuando se usa la DFT, ya que solo estamos tratando con un número finito de muestras. Desafortunadamente, la DFT de una onda cuadrada es una función de tipo sinc (sin (x) / x).

El problema de tener transiciones bruscas en su filtro (salto rápido de 0 a 1 en una muestra) es que esto tiene una respuesta muy larga en el dominio del tiempo, que se está truncando por una ventana cuadrada. Entonces, 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 de ventanas agregando la línea:

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

Después de tomar el IFFT, obtenemos esta respuesta:

respuesta de frecuencia http://www.freeimagehosting.net/uploads/944da9dd93.png

Así que recomendaría intentar implementar el método de diseño de ventanas ya que es bastante simple (hay mejores formas, pero se complican más). Ya que está implementando un ecualizador, asumo que desea poder cambiar las atenuaciones sobre la marcha, por lo que sugeriría calcular y almacenar el filtro en el dominio de la frecuencia siempre que haya un cambio en los parámetros, y luego simplemente puede aplicarlo a cada búfer de audio de entrada tomando los pies por pulgada del búfer de entrada, multiplicando por sus muestras de filtro de dominio de frecuencia y luego ejecutando el ifft para volver al dominio de tiempo. Esto será mucho más eficiente que todas las ramificaciones que realice para cada muestra.

Answer 2

Hay dos problemas: la forma en que usa 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 usa la Transformada Discreta de Fourier (DFT) (implementada con un algoritmo de Transformada Rápida de Fourier para la velocidad), en realidad calcula una versión muestreada del verdadero espectro. 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 la sintaxis de Matlab / Octave:

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:

La falla al principio del bloque no es lo que esperamos en absoluto. Pero si consideras fft(x) , tiene sentido. La transformada discreta de Fourier asume que la señal es periódica dentro del bloque de transformación. Hasta donde sabe el DFT, pedimos la transformación de un período de esto:

Esto lleva a la primera consideración importante cuando se filtra con DFT: en realidad está implementando convolución circular , no convolución lineal. Por lo tanto, la "falla" en el primer gráfico no es realmente una falla cuando se consideran las matemáticas. Entonces, la pregunta se convierte en: ¿hay una manera de evitar la periodicidad? La respuesta es sí: usar el procesamiento de superposición-guardar . Esencialmente, usted calcula los productos N-largos como se indicó anteriormente, 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 :

Esta imagen tiene sentido: esperamos un transitorio de inicio del filtro, luego el resultado se asienta en la respuesta sinusoidal de estado estable. 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 tema: el filtro en particular. En su bucle, crea un filtro cuyo espectro es esencialmente H = [0 (1:511)/512 1 (511:-1:1)/512]''; Si lo haces hraw = real(ifft(H)); plot(hraw) hraw = real(ifft(H)); plot(hraw) , se obtiene:

Es difícil de ver, pero hay un montón de puntos que no son cero en el extremo izquierdo del gráfico, y luego un montón más en el extremo derecho. Usando la función freqz incorporada de freqz para ver la respuesta de frecuencia que vemos (haciendo freqz(hraw) ):

La respuesta de magnitud tiene muchas ondulaciones desde la envolvente del paso alto hasta cero. Una vez más, la periodicidad inherente a la DFT está en funcionamiento. En lo que respecta a la DFT, hraw repite una y otra vez. Pero si tomas un período de hraw , como freqz hace freqz , su espectro es bastante diferente del de las versiones periódicas.

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

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

Tenga en cuenta que, por lo general, antes de calcular el DFT de la almohadilla h , lo dejé solo en el bucle para compararlo más fácilmente con el original.

Answer 3

Primero, sobre la normalización: eso es un problema conocido (no). La DFT / IDFT requeriría un factor 1 / sqrt (N) (aparte de los factores coseno / seno estándar) en cada uno (directo e inverso) para que sean simétricos y realmente 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 esté al tanto y normalice según lo prefiera. See

Segundo: en una (digamos) DFT de 16 puntos, lo que usted llama la celda 0 correspondería a la frecuencia cero (DC), la frecuencia baja 1 de bin ... la frecuencia media de bin 4 , la bandeja 8 a la frecuencia más alta y las ubicaciones 9. .15 a las "frecuencias negativas". En tu ejemplo, entonces, el contenedor 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 en las frecuencias bajas" . ¿Cómo observas eso?

Answer 4

Su principal problema es que las frecuencias no están bien definidas en intervalos de tiempo cortos . Esto es particularmente cierto para las frecuencias bajas, por lo que es más notorio el problema.

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

Por ejemplo, tomando algunos números razonables: comienzo con una forma de onda a 27,5 Hz (A0 en un piano), digitalizada a 44100 Hz, se verá así (donde la parte roja tiene una longitud de 1024 muestras):

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

Así que primero comenzaremos con un paso bajo de 40Hz. Entonces, como 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? Mal, mal, mal , y esto 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 la DFT.

El hecho de que 27.5 Hz no es particularmente significativo en el segmento corto se puede ver al observar el DFT en la siguiente figura. Tenga en cuenta que aunque la DFT (puntos negros) del segmento más largo muestra un pico a 27,5 Hz, el corto (puntos rojos) no lo hace.

texto alternativo http://i50.tinypic.com/14w6luw.png

Claramente, luego el filtrado por debajo de 40Hz, solo capturará el desplazamiento de CC, y el resultado del filtro de paso bajo de 40Hz se muestra en verde a continuación.

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

La curva azul (tomada con un corte de 200 Hz) está comenzando a coincidir mucho mejor. Pero tenga en cuenta que no son las bajas frecuencias las que hacen que coincida bien, sino la inclusión de altas frecuencias. No es hasta que incluimos todas las frecuencias posibles en el segmento corto, hasta 22KHz, 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 no es una onda sinusoidal de 27,5 Hz, y su DFT no tiene mucho que ver con 27,5 Hz.