DSP - filtrage dans le domaine de fréquence via FFT

Question

Je joue autour d'un peu avec la mise en œuvre Exocortex de la FFT, mais je vais avoir quelques problèmes.

Chaque fois que je modifie les amplitudes des bandes de fréquence avant d'appeler le iFFT le signal résultant contient des bruits parasites, en particulier lorsque les fréquences basses sont présentes dans le signal (comme des tambours ou basses). Cependant, cela ne se produit pas si j'Atténuer tous les bacs par le même facteur.

Permettez-moi de mettre un exemple de la mémoire tampon de sortie d'une FFT 4-échantillon:

// 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 sortie se compose de paires de flotteurs, chacune représentant des parties réelles et imaginay d'un seul bac. Donc, bin 0 (index de tableau 0, 1) représenterait les parties réelles et imaginaires de la fréquence DC. Comme vous pouvez le voir, les bacs 1 et 3 ont tous deux les mêmes valeurs, (sauf pour le signe de la partie Im), donc je suppose que bin 3 est la première fréquence négative, et enfin des index (4, 5) serait le dernier positif bin fréquence.

Ensuite, pour atténuer le bac de fréquence 1 ce que je fais:

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

Pour les essais réels que je utilise une FFT de 1024 longueur et je toujours fournir tous les échantillons donc pas 0-remplissage est nécessaire.

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

Il est évident que je fais quelque chose de mal mais ne peut pas comprendre ce que.

Je ne veux pas utiliser la sortie FFT comme un moyen pour générer un ensemble de coefficients FIR depuis que je suis en train de mettre en œuvre un égaliseur dynamique très basique.

Quelle est la bonne façon de filtre dans le domaine de fréquence? ce que je suis absent?

En outre, est-il vraiment nécessaire d'atténuer les fréquences négatives aussi bien? Je l'ai vu une implémentation FFT où neg. les valeurs de fréquence sont mis à zéro avant la synthèse.

Merci à l'avance.

Answer 1

Il y a deux problèmes: la façon dont vous utilisez la FFT, et le filtre particulier

.

Le filtrage est traditionnellement mis en œuvre comme convolution dans le domaine temporel. Vous avez raison que la multiplication des spectres des signaux d'entrée et de filtrage est équivalent. Toutefois, lorsque vous utilisez la transformée de Fourier discrète (DFT) (mis en œuvre avec un algorithme de transformation de Fourier rapide pour la vitesse), on calcule en fait une version échantillonnée du spectre vrai. Cela a beaucoup de conséquences, mais le plus pertinent pour le filtrage est l'implication que le signal de domaine temporel est périodique.

Voici un exemple. Considérons un signal d'entrée sinusoïdal x avec 1,5 cycles de la période, et un h simple filtre passe-bas. Dans la syntaxe 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);

Et voici le graphique:

Le pépin au début du bloc est pas ce que nous attendons du tout. Mais si l'on considère fft(x), il est logique. La transformée de Fourier discrète suppose que le signal est périodique dans le bloc de transformée. En ce qui concerne le DFT sait, nous avons demandé la transformation d'une période de celle-ci:

Cela conduit à la première considération importante lors du filtrage avec DFT: vous implémentez fait circulaire convolution , non convolution linéaire. Ainsi, le « pépin » dans le premier graphique n'est pas vraiment un petit problème quand on considère les mathématiques. Donc, la question devient alors: est-il un moyen de contourner la périodicité? La réponse est oui: utiliser chevauchement sauvegarde traitement . Essentiellement, vous calculez les produits N-longs comme ci-dessus, mais seulement garder N / 2 points.

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

Et voici le graphique de ycorrect:

Cette image est logique - nous nous attendons à un transitoire de démarrage du filtre, le résultat se installe dans la réponse sinusoïdale de l'état d'équilibre. Notez que maintenant x peut être arbitrairement longue. La limitation est Nproc > 2*min(length(x),length(h)).

sur la deuxième question: le filtre particulier. Dans votre boucle, vous créez un filtre qui est le spectre est essentiellement H = [0 (1:511)/512 1 (511:-1:1)/512]'; Si vous hraw = real(ifft(H)); plot(hraw), vous obtenez:

Il est difficile de voir, mais il y a un tas de points non nuls au bord extrême gauche du graphique, puis un groupe plus au bord extrême droite. L'utilisation d'Octave fonction freqz intégré pour regarder la réponse en fréquence, nous voyons (en faisant freqz(hraw)):

La réponse de magnitude a beaucoup d'ondulations de l'enveloppe passe-haut à zéro. Encore une fois, la périodicité inhérente à la DFT est au travail. En ce qui concerne la TFD concerne, répète hraw encore et encore. Mais si vous prenez une période de hraw, comme freqz fait, son spectre est tout à fait différent de la version de périodique.

Alors, nous allons définir un nouveau signal: Nous hrot = [hraw(513:end) ; hraw(1:512)]; simplement la sortie alternons brute DFT pour le rendre continu à l'intérieur du bloc. Maintenant, regardons la réponse en fréquence en utilisant freqz(hrot):

Beaucoup mieux. L'enveloppe souhaitée est là, sans que toutes les ondulations. Bien sûr, la mise en œuvre est pas si simple maintenant, vous devez faire une multiplication complète complexe par fft(hrot) plutôt que chaque mise à l'échelle bin complexe, mais au moins vous aurez la bonne réponse.

Notez que pour la vitesse, vous souhaitez généralement pré-calculer la DFT du h rembourré, je l'ai laissé seul dans la boucle à plus facily comparer avec l'original.

Answer 2

Votre principal problème est que fréquences ne sont pas bien définies sur des intervalles de temps courts . Cela est particulièrement vrai pour les basses fréquences, ce qui est la raison pour laquelle vous remarquez le problème le plus là-bas.

Par conséquent, lorsque vous prenez des segments très courts du train sonore, puis vous filtrez ceux-ci, les segments filtrés filtre wont d'une manière qui produit une forme d'onde continue, et vous entendez les sauts entre les segments et c'est ce qui génère le clics ici.

Par exemple, en prenant quelques chiffres raisonnables: Je commence par une forme d'onde à 27,5 Hz (A0 sur un piano), numérisé à 44100 Hz, il ressemblera à ceci (où la partie rouge est longue de 1024 échantillons):

Alors d'abord, nous allons commencer par une passe de bas de 40Hz. Donc, puisque la fréquence d'origine est inférieure à 40 Hz, un filtre passe-bas avec une coupure de 40Hz devrait pas vraiment d'effet, et nous allons obtenir une sortie qui correspond presque exactement à l'entrée. Droite? Faux, faux, faux - et cela est essentiellement le noyau de votre problème. Le problème est que pour les courtes sections idée de 27,5 Hz est pas clairement définie, et ne peut pas être bien représentés dans la DFT.

Ce 27,5 Hz n'est pas particulièrement significatif dans le court segment peut être vu en regardant la TFD dans la figure ci-dessous. Notez que bien que le DFT de segment plus long (points noirs) montre un pic à 27,5 Hz, le court (points rouges) ne fonctionne pas.

De toute évidence, puis filtration en dessous de 40 Hz, va simplement saisir le décalage en courant continu, et le résultat du filtre passe-bas 40Hz est représenté en vert ci-dessous.

La courbe bleue (prise avec un seuil de coupure de 200 Hz) est en train de faire correspondre beaucoup mieux. Mais notez que ce ne sont pas les basses fréquences qui font correspondre à bien résisté, mais l'inclusion des hautes fréquences. Ce n'est pas jusqu'à ce que nous incluons toutes les fréquences possibles dans le court segment, jusqu'à 22KHz que nous obtenons finalement une bonne représentation de l'onde sinusoïdale d'origine.

La raison de tout cela est qu'un petit segment d'une onde sinusoïdale 27,5 Hz est pas une onde sinusoïdale 27,5 Hz, et il DFT n'a pas grand-chose à voir avec 27,5 Hz.

Answer 3

Êtes-vous la valeur de atténue l'échantillon de fréquence DC à zéro? Il semble que vous n'êtes pas du tout atténuante dans votre exemple. Puisque vous implémentez un filtre passe-haut, vous devez définir la valeur DC à zéro ainsi.

Ceci expliquerait la distorsion basse fréquence. Vous avez beaucoup d'ondulation dans la réponse en fréquence aux basses fréquences si cette valeur DC est non nulle en raison de la grande transition.

est un exemple ici dans Matlab / Octave pour montrer ce qui pourrait se passer:

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

Notez que dans mon code, je crée un exemple de la valeur DC étant non nulle, suivie d'un changement brusque à zéro, puis une montée en puissance. Je prends alors la IFFT pour transformer dans le domaine temporel. Puis-je effectuer une fft zéro rembourré (qui se fait automatiquement par Matlab lorsque vous passez dans une taille de fft plus grand que le signal d'entrée) sur ce signal dans le domaine temporel. Le remplissage avec des zéros dans les résultats dans le domaine temporel en interpolation dans le domaine des fréquences. Avec cela, nous pouvons voir comment le filtre réagira entre les échantillons de filtre.

L'une des choses les plus importantes à retenir est que même si vous définissez des valeurs de réponse de filtre à des fréquences données en atténuant les sorties du DFT, ce ne garantit rien pour les fréquences qui se produisent entre les points d'échantillonnage. Cela signifie que la plus abrupte vos changements, plus dépassement et oscillation entre les échantillons se produiront.

Maintenant, pour répondre à votre question sur la façon dont devrait être fait ce filtrage. Il y a un certain nombre de façons, mais l'un des plus faciles à mettre en œuvre et à comprendre est la méthode de conception de fenêtre. Le problème avec votre conception actuelle est que la largeur de transition est énorme. La plupart du temps, vous voulez aussi rapide des transitions que possible, avec aussi peu d'entraînement possible.

Dans le code suivant, je vais créer un filtre idéal et afficher la réponse:

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

Notez qu'il ya beaucoup d'oscillation provoquée par les changements brusques.

Le FFT ou transformée de Fourier discrète est une version échantillonnée de la transformée de Fourier. La transformée de Fourier est appliquée à un signal sur la plage continue de -infinity à l'infini tandis que le DFT est appliqué sur un nombre fini d'échantillons. Cela entraîne à effet dans un fenêtrage carré (troncature) dans le domaine temporel lors de l'utilisation du DFT puisque nous ne traitons un nombre fini d'échantillons. Malheureusement, la DFT d'une onde carrée est une fonction de type sinus cardinal (sin (x) / x).

Le problème d'avoir des transitions nettes dans le filtre (saut rapide de 0 à 1 dans un seul échantillon) est que ceci a une réponse très longtemps dans le domaine temporel, qui est tronquée par une fenêtre carrée. Donc, pour aider à minimiser ce problème, on peut multiplier le signal dans le domaine temporel par une fenêtre plus progressive. Si nous multiplions une fenêtre de Hanning en ajoutant la ligne:

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

après avoir pris le IFFT, nous obtenons cette réponse:

Je recommande donc d'essayer de mettre en œuvre la méthode de conception de fenêtre, car il est assez simple (il y a de meilleures façons, mais ils deviennent plus compliquées). Puisque vous implémentez un égaliseur, je suppose que vous voulez être en mesure de changer les atténuations à la volée, je vous conseille donc le calcul et le stockage du filtre dans le domaine de fréquence chaque fois qu'il ya un changement dans les paramètres, et vous pouvez simplement l'appliquer mettre en mémoire tampon chaque audio d'entrée en prenant la FFT du tampon d'entrée, en multipliant par les échantillons de filtre de domaine de fréquence, et en effectuant ensuite laIFFT pour revenir au domaine temporel. Ce sera beaucoup plus efficace que tous les vous faites de branchement pour chaque échantillon.

Answer 4

Tout d'abord, au sujet de la normalisation: qui est un problème connu (non). Le DFT / IDFT nécessiterait un facteur 1 / sqrt (N) pour les rendre simmetric et vraiment inversible (en dehors de la norme / cosinus facteurs sinus) dans chacun (direct l'inverse). Une autre possibilité consiste à diviser l'un d'entre eux (le direct ou l'inverse) par N , ceci est une question de commodité et le goût. Souvent, les routines FFT ne fonctionnent pas cette normalisation, l'utilisateur devrait être au courant et normalise comme il préfère. Voir

Deuxièmement: dans un (disons) 16 points DFT, ce que vous appelez bin 0 correspondrait à la fréquence zéro (DC), bin 1 FREQ .. . bin 4 freq moyenne, bin 8 à la plus haute fréquence et bacs 9 ... 15 pour les "fréquences négatives". Vous par exemple, puis, bin 1 est en fait à la fois la fréquence basse et moyenne fréquence. En dehors de cette considération, il n'y a rien sur le plan conceptuel mal dans votre « péréquation ». Je ne comprends pas ce que vous entendez par « le signal est déformé aux basses fréquences » . Comment observez-vous?