DSP - Filtraggio nel dominio della frequenza tramite FFT

Question

Ho giocato in giro un po 'con l'attuazione Exocortex della FFT, ma sto avendo alcuni problemi.

Ogni volta a modificare le ampiezze dei bin di frequenza prima di chiamare l'IFFT il segnale risultante contiene alcuni click e pop, specialmente quando le basse frequenze sono presenti nel segnale (come tamburi o bassi). Tuttavia, questo non succede se io attenuare tutti i bidoni per lo stesso fattore.

Mettiamola un esempio del buffer di uscita di un FFT 4-campione:

// 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'uscita è composto da coppie di galleggianti, ognuno dei quali rappresenta la parte reale e imaginay di un unico raccoglitore. Quindi, bin 0 (indici di matrice 0, 1) rappresenterebbe la parte reale e immaginaria della frequenza DC. Come si può vedere, bidoni 1 e 3 hanno entrambi gli stessi valori, (tranne per il segno della parte Im), quindi credo che bin 3 è la prima frequenza negativo, e, infine, gli indici (4, 5) sarebbe stata l'ultima positiva frequenza bin.

Poi, per attenuare il bidone frequenza di 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 le prove effettive sto usando una FFT a 1024 di lunghezza e ho sempre fornire tutti i campioni in modo non è necessaria alcuna 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 utilizzare l'uscita FFT come un mezzo per generare un insieme di coefficienti FIR visto che sto cercando di attuare un equalizzatore dinamica molto di base.

Qual è il modo corretto di filtro nel dominio della frequenza? quello che mi manca?

Inoltre, è davvero necessaria per attenuare le frequenze negative pure? Ho visto un'implementazione FFT dove neg. valori di frequenza sono azzerati prima della sintesi.

Grazie in anticipo.

Answer 1

Ci sono due problemi: il modo in cui si utilizza la FFT, e il filtro particolare

.

Il filtraggio è tradizionalmente implementato come convoluzione nel dominio del tempo. Hai ragione che moltiplicando gli spettri dei segnali di ingresso e di filtro è equivalente. Tuttavia, quando si utilizza la Discrete Fourier Transform (DFT) (implementato con una Fast Fourier Transform algoritmo per la velocità), è in realtà calcolare una versione campionata del vero spettro. Questo ha un sacco di implicazioni, ma il più rilevante per il filtraggio è l'implicazione che il segnale nel dominio del tempo è periodico.

Ecco un esempio. Si consideri un x segnale di ingresso sinusoidale con 1,5 cicli nel periodo, e una semplice h filtro passa basso. Nella sintassi 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:

Il problema tecnico all'inizio del blocco non è quello che ci aspettiamo affatto. Ma se si considera fft(x), ha senso. Discrete Fourier Transform assume il segnale è periodico all'interno del blocco di trasformazione. Per quanto riguarda il DFT sa, abbiamo chiesto la trasformata di un periodo di questo:

Questo porta alla prima considerazione importante quando si filtra con DFTS: in realtà si implementa circolare convoluzione , convoluzione non lineare. Così il "problema tecnico" nel primo grafico non è davvero un problema tecnico se si considera la matematica. Allora la domanda diventa: c'è un modo per aggirare la periodicità? La risposta è sì: l'uso sovrapposizione-Save elaborazione . In sostanza, si calcola prodotti N-lunghi come sopra, ma solo tenere 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:

Questa immagine ha un senso - ci aspettiamo un transitorio di avvio dal filtro, poi le deposita risultato nella risposta sinusoidale a regime. Si noti che ora x può essere arbitrariamente lungo. La limitazione è Nproc > 2*min(length(x),length(h)).

Ora sul secondo problema: il filtro particolare. Nel suo ciclo, si crea un filtro che è lo spettro è essenzialmente H = [0 (1:511)/512 1 (511:-1:1)/512]'; Se lo fai hraw = real(ifft(H)); plot(hraw), si ottiene:

È difficile vedere, ma ci sono un gruppo di non-zero punti al bordo sinistro del grafico, e poi un gruppo di più al bordo destro. Usando di Octave funzione integrata freqz a guardare la risposta in frequenza che vediamo (facendo freqz(hraw)):

La risposta in ampiezza ha un sacco di increspature dal passa-alto busta fino a zero. Anche in questo caso, la periodicità inerente alla DFT è al lavoro. Per quanto riguarda il DFT è interessato, ripete hraw più e più volte. Ma se si prende un periodo di hraw, come freqz fa, il suo spettro è molto diversa da quella della versione periodica.

Così definiamo un nuovo segnale: hrot = [hraw(513:end) ; hraw(1:512)]; Abbiamo semplicemente ruotare l'output grezzo DFT per renderlo continuo all'interno del blocco. Ora diamo un'occhiata alla risposta di frequenza utilizzando freqz(hrot):

Molto meglio. La busta desiderato è lì, senza tutte le increspature. Naturalmente, l'implementazione non è così semplice ora, si ha a che fare un complesso completo moltiplicare per fft(hrot) piuttosto che scalare ogni bin complesso, ma almeno si otterrà la risposta giusta.

Si noti che per la velocità, si sarebbe di solito pre-calcolare la DFT del h imbottita, ho lasciato da solo nel ciclo a più easily confrontare con l'originale.

Answer 2

Il problema principale è che i le frequenze non sono ben definiti su intervalli di tempo brevi . Ciò è particolarmente vero per le basse frequenze, che è il motivo per cui si nota il problema più lì.

Di conseguenza, quando si prende segmenti davvero brevi dal treno suono, e quindi di filtrare questi, i segmenti filtrati filtro di solito in un modo che produce una forma d'onda continua, e si sente il salti tra i segmenti e questo è ciò che genera la click qui.

Ad esempio, prendendo alcuni numeri ragionevoli: Comincio con una forma d'onda a 27,5 Hz (A0 su un pianoforte), digitalizzato a 44100 Hz, sarà simile a questa (dove la parte rossa è lungo 1024 campioni):

Quindi, prima inizieremo con un passa-basso di 40Hz. Quindi dal momento che la frequenza originale è inferiore a 40Hz, un filtro passa-basso con una 40Hz cut-off dovrebbe in realtà non ha alcun effetto, e noi ottenere un output che corrisponde quasi esattamente l'input. Giusto? Sbagliato, sbagliato, sbagliato - e questo è fondamentalmente il nocciolo del problema. Il problema è che per brevi sezioni idea di 27,5 Hz non è chiaramente definito, e non può essere rappresentato bene nel DFT.

Quello 27,5 Hz non è particolarmente significativa nel breve segmento può essere visto, cercando in DFT nella figura sottostante. Nota che, anche se DFT del segmento più lungo (punti neri) mostra un picco a 27,5 Hz, quello corto (punti rossi) non lo fa.

Chiaramente, quindi filtrando sotto 40Hz, sarà solo catturare l'offset DC, e il risultato del filtro passa-basso 40 Hz è mostrato in verde qui sotto.

La curva blu (prese con un 200 Hz cut-off) sta iniziando a corrispondere meglio. Ma nota che non è le basse frequenze che stanno rendendo all'altezza bene, ma l'inserimento di alte frequenze. Non è fino a quando includiamo tutte le frequenze possibili nel breve segmento, fino a 22KHz che finalmente otteniamo una rappresentazione buona della sinusoide originale.

La ragione di tutto questo è che un piccolo segmento di un'onda sinusoidale 27,5 Hz è non un'onda sinusoidale 27,5 Hz, e di DFT non ha molto a che fare con 27,5 Hz.

Answer 3

Stai interferisca con il valore del campione di frequenza DC a zero? Sembra che non si sta attenuando affatto nel tuo esempio. Poiché si implementa un filtro passa alto, è necessario impostare il valore DC a zero pure.

Questo spiegherebbe la distorsione delle basse frequenze. Si avrebbe un sacco di ripple nella risposta in frequenza alle basse frequenze se tale valore DC è diverso da zero a causa della grande transizione.

Ecco un esempio in MATLAB / Octave per dimostrare ciò che 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');

Si noti che nel mio codice, sto creando un esempio del valore DC essendo diverso da zero, seguito da un brusco cambiamento a zero, e poi una rampa in salita. Ho poi prendo l'IFFT di trasformare nel dominio del tempo. Poi eseguo un FFT zeri (che viene fatto automaticamente dal MATLAB quando si passa a una dimensione più grande di FFT il segnale in ingresso) su quel segnale nel dominio del tempo. Il riempimento di zeri nei risultati nel dominio del tempo in interpolazione nel dominio della frequenza. Usando questo, possiamo vedere come il filtro risponderà tra i campioni di filtro.

Una delle cose più importanti da ricordare è che anche se si sta impostando i valori di risposta del filtro a date frequenze attenuando le uscite del DFT, questo garantisce nulla per frequenze intermedie comprese punti di campionamento. Questo significa che il più brusco le modifiche, si verificheranno più overshoot e oscillazione tra i campioni.

Ora per rispondere alla tua domanda su come questo filtraggio dovrebbe essere fatto. Ci sono un certo numero di modi, ma uno dei più facili da implementare e capire è il metodo della finestra di progettazione. Il problema con il vostro disegno attuale è che la larghezza di transizione è enorme. Il più delle volte, si vuole il più velocemente di transizioni il più possibile, con il minor ondulazione possibile.

Nel codice seguente, ho creerà un filtro ideale e visualizzare 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'); 

Si noti che c'è un sacco di oscillazione causata dai cambiamenti bruschi.

La FFT o Discrete Fourier Transform è una versione campionata della trasformata di Fourier. La trasformata di Fourier è applicato ad un segnale nell'intervallo -infinito continuo all'infinito mentre la DFT viene applicato sopra un numero finito di campioni. Questo a effetto in un finestre quadrato (troncamento) nel dominio del tempo quando si utilizza il DFT dal momento che si tratta solo di un numero finito di campioni. Purtroppo, la DFT di un'onda quadra è una funzione sinc tipo (sin (x) / x).

Il problema di avere transizioni taglienti nel filtro (salto rapido 0-1 in un campione) è che questo ha una lunga risposta nel dominio del tempo, che viene troncato da una finestra quadrata. Quindi, per ridurre al minimo questo problema, siamo in grado di moltiplicare il segnale nel dominio del tempo da una finestra più graduale. Se moltiplichiamo una finestra Hanning aggiungendo la riga:

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

dopo aver preso l'IFFT, otteniamo questa risposta:

Quindi mi sento di raccomandare cercando di attuare il metodo della finestra di progettazione dal momento che è abbastanza semplice (ci sono modi migliori, ma ottengono più complicato). Dal momento che si sta implementando un equalizzatore, presumo si vuole essere in grado di cambiare le attenuazioni al volo, quindi vorrei suggerire calcolare e memorizzare il filtro nel dominio della frequenza ogni volta che c'è un cambiamento di parametri, e allora si può solo applicare ad ogni ingresso audio tamponare prendendo la FFT del buffer di ingresso, moltiplicando per i vostri campioni di filtro nel dominio della frequenza, e quindi eseguire laIFFT Per tornare al dominio del tempo. Questo sarà molto più efficiente di tutte le ramificazioni che si sta facendo per ogni campione.

Answer 4

In primo luogo, circa la normalizzazione: che è un noto (non) problema. La DFT / IDFT richiederebbe un fattore di 1 / sqrt (N) (a parte il coseno di serie / fattori sinusoidali) in ognuno (dirigere un inversa) per renderli simmetrica e veramente invertibile. Un'altra possibilità è quella di dividere uno di loro (la diretta o l'inverso) da N , questa è una questione di convenienza e gusto. Spesso le routine FFT non eseguono questa normalizzazione, l'utente dovrebbe essere a conoscenza di esso e normalizzare come preferisce. Vedi

In secondo luogo: in un (diciamo) 16 punti DFT, quello che si chiama il bin 0 corrisponderebbe alla frequenza zero (DC), bin 1 bassa frequenza .. . bin 4 di media freq, bin 8 per la più alta frequenza e bidoni 9 ... 15 per le "frequenze negative". In te esempio, allora, bin 1 è in realtà sia la frequenza bassa e media frequenza. A parte questa considerazione, non c'è nulla di concettualmente sbagliato nel vostro "perequazione". Non capisco cosa intendi per "il segnale viene distorta a basse frequenze" . Come si osserva che?