DSP: filtraggio nel dominio della frequenza tramite FFT

Ho giocato un po ‘con l’implementazione Exocortex della FFT, ma ho qualche problema.

Ogni volta che modifico le ampiezze dei raccoglitori di frequenza prima di chiamare l’iFFT, il segnale risultante contiene alcuni clic e pop, specialmente quando nel segnale sono presenti basse frequenze (come batteria o bassi). Tuttavia, ciò non accade se attenuo tutti i contenitori per lo stesso fattore.

Consentitemi di mettere un esempio del buffer di output di un FFT a 4 campioni:

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

L’output è composto da coppie di float, ognuno dei quali rappresenta le parti reali e immaginarie di un singolo bin. Quindi, bin 0 (array indexes 0, 1) rappresenterebbe le parti reali e immaginarie della frequenza DC. Come puoi vedere, i bin 1 e 3 hanno entrambi gli stessi valori, (eccetto il segno della parte Im), quindi suppongo che bin 3 sia la prima frequenza negativa, e infine gli indici (4, 5) sarebbero gli ultimi positivi scomparto di frequenza.

Quindi per attenuare il bin di frequenza 1 questo è quello che faccio:

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

Per i test effettivi sto usando un FFT da 1024 lunghezze e fornisco sempre tutti i campioni, quindi non è necessario uno 0-padding.

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

Ovviamente sto facendo qualcosa di sbagliato ma non riesco a capire cosa.

Non voglio usare l’output FFT come mezzo per generare un set di coefficienti FIR poiché sto cercando di implementare un equalizzatore dinamico molto semplice.

Qual è il modo corretto di filtrare nel dominio della frequenza? cosa mi manca?

Inoltre, è davvero necessario attenuare anche le frequenze negative? Ho visto un’implementazione FFT dove neg. i valori di frequenza vengono azzerati prima della sintesi.

Grazie in anticipo.

Ci sono due problemi: il modo in cui usi la FFT e il filtro particolare.

Il filtraggio è tradizionalmente implementato come convoluzione nel dominio del tempo. Hai ragione che moltiplicando gli spettri dell’input e i segnali del filtro è equivalente. Tuttavia, quando si utilizza la Discrete Fourier Transform (DFT) (implementata con un algoritmo Fast Fourier Transform per la velocità), si calcola effettivamente una versione campionata dello spettro reale. Ciò ha molte implicazioni, ma quella più rilevante per il filtraggio è l’implicazione che il segnale del dominio del tempo sia periodico.

Ecco un esempio. Considerare un segnale di ingresso sinusoidale x con 1,5 cicli nel periodo e un filtro passa-basso semplice h . Nella syntax 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); 

Ed ecco il grafico: IFFT del prodotto

Il problema tecnico all’inizio del blocco non è quello che ci aspettiamo. Ma se consideri fft(x) , ha senso. La Trasformata di Fourier discreta presuppone che il segnale sia periodico all’interno del blocco di trasformazione. Per quanto ne sa la DFT, abbiamo chiesto la trasformazione di un periodo di questo: Input aperiodico a DFT

Questo porta alla prima considerazione importante quando si filtra con DFT: si sta effettivamente implementando la convoluzione circolare , non la convoluzione lineare. Quindi il “glitch” nel primo grafico non è proprio un problema tecnico quando si considera la matematica. Quindi la domanda diventa: c’è un modo per aggirare la periodicità? La risposta è sì: usa l’ elaborazione di sovrapposizione-salvataggio . In sostanza, si calcolano prodotti N-long come sopra, ma si mantengono solo N / 2 punti.

 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 

Ed ecco il grafico di ycorrect : ycorrect

Questa immagine ha un senso – ci aspettiamo un transitorio di avvio dal filtro, quindi il risultato si assesta nella risposta sinusoidale allo stato stazionario. Nota che ora x può essere arbitrariamente lungo. La limitazione è Nproc > 2*min(length(x),length(h)) .

Ora sul secondo problema: il filtro particolare. Nel tuo ciclo, crei un filtro il cui spettro è essenzialmente H = [0 (1:511)/512 1 (511:-1:1)/512]'; Se fai hraw = real(ifft(H)); plot(hraw) hraw = real(ifft(H)); plot(hraw) , ottieni: hraw

È difficile da vedere, ma ci sono una serie di punti diversi da zero all’estremità sinistra del grafico, e poi un mucchio in più all’estremità destra. Usando la funzione freqz integrata di freqz per osservare la risposta in frequenza che vediamo (facendo freqz(hraw) ): freqz (hraw)

La risposta di grandezza ha un sacco di increspature dalla busta passa-alto fino a zero. Ancora una volta, la periodicità inerente alla DFT è al lavoro. Per quanto riguarda la DFT, l’ hraw ripete all’infinito. Ma se prendi un periodo di hraw , come fa freqz , il suo spettro è molto diverso da quello della versione periodica.

Quindi definiamo un nuovo segnale: hrot = [hraw(513:end) ; hraw(1:512)]; hrot = [hraw(513:end) ; hraw(1:512)]; Semplicemente ruotiamo l’output DFT grezzo per renderlo continuo all’interno del blocco. Ora diamo un’occhiata alla risposta in frequenza usando freqz(hrot) : freqz (hrot)

Molto meglio. La busta desiderata è lì, senza tutte le increspature. Naturalmente, l’implementazione non è così semplice ora, devi fare un multiplo complesso completo per fft(hrot) invece di ridimensionare ogni bin complesso, ma almeno avrai la risposta giusta.

Si noti che per la velocità, di solito si pre-calcolare la DFT della h imbottita, l’ho lasciato da solo nel ciclo per confrontare più facilmente con l’originale.

Il tuo problema principale è che le frequenze non sono ben definite su intervalli di tempo brevi . Questo è particolarmente vero per le basse frequenze, motivo per cui notate il problema più spesso.

Pertanto, quando si estrae segmenti molto brevi dal treno sonoro, e quindi si filtrano questi, i segmenti filtrati non filtrano in un modo che produce una forma d’onda continua, e si sentono i salti tra i segmenti e questo è ciò che genera i clic qui .

Ad esempio, prendendo alcuni numeri ragionevoli: inizio con una forma d’onda a 27,5 Hz (A0 su un piano), digitalizzata a 44100 Hz, sarà simile a questa (dove la parte rossa è lunga 1024 campioni):

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

Quindi per prima cosa inizieremo con un passa-basso di 40Hz. Quindi, poiché la frequenza originale è inferiore a 40 Hz, un filtro passa-basso con un cut-off a 40 Hz non dovrebbe avere alcun effetto, e otterremo un risultato che corrisponde esattamente all’ingresso. Destra? Sbagliato, sbagliato, sbagliato – e questo è fondamentalmente il nocciolo del tuo problema. Il problema è che per le brevi sezioni l’ idea di 27,5 Hz non è chiaramente definita e non può essere rappresentata bene nella DFT.

Quella 27.5 Hz non è particolarmente significativa nel segmento corto può essere visto guardando la DFT nella figura sottostante. Si noti che sebbene il DFT (punti neri) del segmento più lungo mostri un picco a 27,5 Hz, quello corto (punti rossi) no.

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

Chiaramente, quindi il filtraggio sotto i 40Hz, catturerà solo l’offset DC, e il risultato del filtro passa-basso a 40Hz è mostrato in verde sotto.

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

La curva blu (presa con un cutoff a 200 Hz) inizia a corrispondere molto meglio. Ma nota che non sono le basse frequenze a farla combaciare bene, ma l’inclusione delle alte frequenze. Non è fino a quando non includiamo tutte le frequenze possibili nel segmento corto, fino a 22 KHz che finalmente otteniamo una buona rappresentazione dell’onda sinusoidale originale.

La ragione di tutto ciò è che un piccolo segmento di un’onda sinusoidale di 27,5 Hz non è un’onda sinusoidale di 27,5 Hz, ed è che DFT non ha molto a che fare con 27,5 Hz.

Stai attenuando il valore del campione di frequenza DC a zero? Sembra che non lo stiate attenuando affatto nel tuo esempio. Poiché stai implementando un filtro passa-alto, devi anche impostare il valore DC su zero.

Ciò spiegherebbe la distorsione a bassa frequenza. Avresti un sacco di ripple nella risposta in frequenza alle basse frequenze se quel valore DC è diverso da zero a causa della grande transizione.

Ecco un esempio in MATLAB / Octave per dimostrare cosa potrebbe accadere:

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

risposta in frequenza http://www.freeimagehosting.net/uploads/e10109e535.png

Si noti che nel mio codice sto creando un esempio in cui il valore DC è diverso da zero, seguito da un brusco cambiamento a zero e quindi da una rampa di accelerazione. Quindi prendo l’IFFT per trasformarlo nel dominio del tempo. Quindi eseguo un fft a zero (che viene eseguito automaticamente da MATLAB quando passi in una dimensione di fft maggiore del segnale di ingresso) su quel segnale del dominio del tempo. Lo zero padding nel dominio del tempo produce un’interpolazione nel dominio della frequenza. Usando questo, possiamo vedere come il filtro risponderà tra i campioni del filtro.

Una delle cose più importanti da ricordare è che anche se si impostano i valori di risposta del filtro a determinate frequenze attenuando le uscite della DFT, ciò non garantisce nulla per le frequenze che si verificano tra i punti campione. Ciò significa che più cambierete bruscamente, più sovraoscillazioni e oscillazioni tra i campioni si verificheranno.

Ora per rispondere alla tua domanda su come dovrebbe essere fatto questo filtraggio. Esistono diversi modi, ma uno dei più facili da implementare e comprendere è il metodo di progettazione delle windows. Il problema con il tuo design attuale è che la larghezza della transizione è enorme. Il più delle volte, vorrai il più veloce ansible delle transizioni, con il minimo ripple ansible.

Nel prossimo codice creerò un filtro ideale e visualizzerò la risposta:

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

risposta in frequenza http://www.freeimagehosting.net/uploads/c86f5f1700.png

Si noti che vi sono molte oscillazioni causate da cambiamenti improvvisi.

FFT o Discrete Fourier Transform è una versione campionata della Trasformata di Fourier. La trasformata di Fourier viene applicata a un segnale nell’intervallo continuo -infinità all’infinito mentre la DFT viene applicata su un numero finito di campioni. Questo in effetti produce una finestra quadrata (troncatura) nel dominio del tempo quando si usa la DFT poiché si tratta solo di un numero finito di campioni. Sfortunatamente, la DFT di un’onda quadra è una funzione di tipo sinc (sin (x) / x).

Il problema con le transizioni brusche nel filtro (salto rapido da 0 a 1 in un campione) è che questa ha una risposta molto lunga nel dominio del tempo, che viene troncata da una finestra quadrata. Quindi, per ridurre al minimo questo problema, possiamo moltiplicare il segnale del dominio del tempo con una finestra più graduale. Se moltiplichiamo una finestra di blocco aggiungendo la linea:

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

dopo aver preso l’IFFT, otteniamo questa risposta:

risposta in frequenza http://www.freeimagehosting.net/uploads/944da9dd93.png

Quindi consiglierei di provare ad implementare il metodo di progettazione delle windows dato che è abbastanza semplice (ci sono modi migliori, ma diventano più complicati). Dato che stai implementando un equalizzatore, presumo che tu voglia essere in grado di modificare le attenuazioni al volo, quindi suggerirei di calcolare e memorizzare il filtro nel dominio della frequenza ogni volta che c’è un cambiamento nei parametri, e quindi puoi semplicemente applicarlo a ciascun buffer audio di input prendendo il fft del buffer di input, moltiplicando per i campioni del filtro del dominio della frequenza e quindi eseguendo l’ifft per tornare al dominio del tempo. Questo sarà molto più efficiente di tutte le ramificazioni che stai facendo per ogni campione.

Innanzitutto, sulla normalizzazione: si tratta di un problema noto (non). Il DFT / IDFT richiederebbe un fattore 1 / sqrt (N) (a parte i coseno / fattori seno standard) in ciascuno (diretto un inverso) per renderli simmetrici e invertibili. Un’altra possibilità è dividere uno di essi (il diretto o l’inverso) di N , questa è una questione di convenienza e gusto. Spesso le routine FFT non eseguono questa normalizzazione, l’utente dovrebbe esserne a conoscenza e normalizzarsi come preferisce. Vedere

Secondo: in un DFT da 16 punti, quello che chiami il bin 0 corrisponde alla frequenza zero (DC), bin 1 low freq … bin 4 media freq, bin 8 alla frequenza più alta e bin 9 .. .15 alle “frequenze negative”. Nell’esempio, quindi, bin 1 è in realtà sia la bassa frequenza che la frequenza media. A prescindere da questa considerazione, non vi è nulla di concettualmente sbagliato nella vostra “equalizzazione”. Non capisco cosa intendi con “il segnale viene distorto a basse frequenze” . Come lo osservi?