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

Question

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

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

Permítanme poner un ejemplo de la memoria intermedia de salida de una FFT de 4-ejemplo:

// 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. Así, bin 0 (matriz índices 0, 1) representaría las partes real e imaginaria de la frecuencia de CC. Como se puede ver, contenedores de 1 y 3 ambos tienen los mismos valores, (excepto el signo de la parte Im), así que supongo bin 3 es la primera frecuencia negativa, y, finalmente, los índices (4, 5) serían los últimos positivo frecuencia bin.

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 proporcionar todas las muestras de modo que no se necesita 0-relleno.

// 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;
}

Es obvio que estoy haciendo algo mal, pero no puedo entender por qué.

No quiero utilizar la salida de la 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 filtro en el dominio de la frecuencia? lo que me falta?

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

Gracias de antemano.

Answer 1

Hay dos cuestiones: el modo de usar la FFT, y el filtro especial

.

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

Este es un ejemplo. Considere un x señal de entrada sinusoidal con 1,5 ciclos en el período, y un simple h filtro de paso bajo. 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:

El problema técnico al inicio del bloque no es lo que esperamos en absoluto. Pero si se tiene en cuenta fft(x), tiene sentido. La Transformada Discreta de Fourier asume que la señal es periódica en el bloque de transformada. 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: se está implementando en realidad circular convolución , convolución no lineal. Por lo que el "fallo" en el primer gráfico no es realmente un problema técnico si tenemos en cuenta los cálculos. Entonces la pregunta es: ¿hay una manera de evitar la periodicidad? La respuesta es sí: el uso solapamiento-ahorra recursos . En esencia, se calcula productos N-largas que el anterior, pero sólo mantener 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, luego se asienta la resultado en la respuesta sinusoidal de estado estacionario. Tenga en cuenta que ahora x puede ser arbitrariamente larga. La limitación es Nproc > 2*min(length(x),length(h)).

Ahora sobre el segundo problema: el filtro particular. En el bucle, se crea un filtro que es esencialmente espectro se 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 montón más en el extremo derecho. El uso de Octave función integrada de freqz mirar la respuesta de frecuencia vemos (haciendo freqz(hraw)):

La respuesta de magnitud tiene una gran cantidad de ondas de la baja dotación de paso alto a cero. Una vez más, la periodicidad inherente a la DFT es en el trabajo. En lo que se refiere a la DFT, repeticiones hraw una y otra vez. Pero si se toma un período de hraw, como freqz hace, su espectro es muy diferente de la versión del periódico.

Así que vamos a definir una nueva señal: hrot = [hraw(513:end) ; hraw(1:512)]; Nos basta con girar la salida cruda DFT para que sea continua dentro del bloque. Ahora vamos a ver la respuesta de frecuencia usando freqz(hrot):

Mucho mejor. El sobre deseado es allí, sin todas las ondulaciones. Por supuesto, la aplicación no es tan simple ahora, usted tiene que hacer un complejo completo se multiplica por fft(hrot) en lugar de escalar cada bandeja compleja, pero al menos obtendrá la respuesta correcta.

Tenga en cuenta que para la velocidad, usted por lo general pre-calcular la DFT de la h acolchado, me dejó solo en el bucle para más easily comparar con el original.

Answer 2

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

Por lo tanto, cuando se toma segmentos muy cortos del tren de sonido, y luego filtrar estos, los segmentos filtrados costumbre de filtro de una manera que produce una forma de onda continua, y se oye los saltos entre los segmentos y esto es lo que genera la usted hace clic aquí.

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

Así que primero vamos a empezar con un paso bajo de 40 Hz. Así que ya que la frecuencia original es inferior a 40Hz, un filtro de paso bajo con un corte de 40 Hz debe en realidad no tiene ningún efecto, y vamos a tener una salida que coincide casi exactamente la entrada. ¿Correcto? Mal, mal, mal - y esto es básicamente el núcleo de su problema. El problema es que para las secciones cortas de la idea de 27,5 Hz no está claramente definido, y no se puede representar bien en la DFT.

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

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

La curva azul (tomado con un 200 Hz de corte) está empezando a coincidir mucho mejor. Pero tenga en cuenta que no es las bajas frecuencias que están haciendo que se acople bien, pero la inclusión de altas frecuencias. No es hasta que incluimos todas las frecuencias posibles en el segmento corto, hasta 22KHz que por fin tenemos una buena representación de la onda sinusoidal originales.

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

Answer 3

¿Usted está atenuando el valor de la muestra de frecuencia DC a cero? Parece que no está en absoluto atenuante en su ejemplo. Dado que va a implementar un filtro de paso alto, es necesario establecer el valor de CC a cero también.

Esto explicaría la distorsión de baja frecuencia. Usted tiene una gran cantidad de ondulación en la respuesta de frecuencia en las frecuencias bajas si ese valor DC no es cero debido a la gran transición.

Este es 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');

Tenga 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 rampa hacia arriba. entonces tomo la IFFT de transformarse en el dominio del tiempo. Entonces puedo realizar una FFT con relleno de ceros (que se realiza automáticamente por MATLAB cuando se pasa en un tamaño de FFT más grande que la señal de entrada) en la que la señal de dominio de tiempo. El relleno con ceros en el dominio del tiempo como resultado de interpolación en el dominio de la frecuencia. Usando esto, podemos ver cómo va a responder el filtro entre muestras de filtro.

Una de las cosas más importantes a tener en cuenta es que a pesar de que va a configurar los valores de respuesta del filtro en las frecuencias dadas por la atenuación de las salidas de la DFT, esto no garantiza de frecuencias situadas entre los puntos de muestreo. Esto significa que el más abrupto sus cambios, se producirá la mayor exceso y la oscilación entre las muestras.

Ahora, para responder a su pregunta de cómo esta filtración se debe hacer. Hay un número de maneras, pero una de las más fáciles de implementar y entender es el método de diseño de la ventana. El problema con el diseño actual es que el ancho de transición es enorme. La mayoría de las veces, tendrá que tan rápido como sea posible de las transiciones, con la menor ondulación como sea posible.

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

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

La FFT o Transformada Discreta de Fourier es una versión muestreada de la transformada de Fourier. La transformada de Fourier se aplica a una señal sobre la -infinito continua rango hasta el infinito, mientras que la DFT se aplica sobre un número finito de muestras. Esto, a efecto de generar un cuadrado de ventanas (truncamiento) en el dominio del tiempo cuando se usa la DFT ya que sólo se trata de un número finito de muestras. Por desgracia, la DFT de una onda cuadrada es una función de tipo sinc (sin (x) / x).

El problema de tener transiciones bruscas en el filtro (salto rápido de 0 a 1 en una muestra) es que esto tiene un muy largo respuesta en el dominio del tiempo, que se está truncado por una ventana cuadrada. Así que para ayudar a minimizar este problema, podemos multiplicar la señal de dominio de 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:

Por lo que recomiendo tratando de poner en práctica el método de diseño de la ventana, ya que es bastante sencillo (hay mejores maneras, pero ellos se ponen más complicadas). Dado que va a implementar un ecualizador, asumo que quiere ser capaz de cambiar las atenuaciones sobre la marcha, así que yo sugeriría calcular y almacenar el filtro en el dominio de la frecuencia cada vez que hay un cambio en los parámetros, y entonces sólo puedo aplicarlo a cada uno de audio de entrada tampón mediante la adopción de la FFT de la memoria intermedia de entrada, multiplicando por sus muestras de filtros de dominio de frecuencia, y luego realizar laifft para volver al dominio del tiempo. Esto va a ser mucho más eficiente que la totalidad de la ramificación que está haciendo para cada muestra.

Answer 4

En primer lugar, sobre la normalización: que es un (no) problema conocido. La DFT / IDFT requeriría un factor 1 / sqrt (N) (aparte del coseno estándar / factores de seno) en cada uno (dirigir una inversa) para hacerlos simmetric y verdaderamente invertible. Otra posibilidad es dividir una de ellas (la directa o la inversa) por N , esto es una cuestión de conveniencia y gusto. A menudo las rutinas FFT no realizan esta normalización, se espera que el usuario sea consciente de ello y normalizar como él prefiere. Ver

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