DSP - Filtrage dans le domaine fréquentiel par la FFT

Il y a deux problèmes: l'utilisation de la FFT, et le filtre particulier.

De filtrage est traditionnellement mis en œuvre en tant que produit de convolution dans le domaine temporel. Vous avez raison que le fait de multiplier les spectres de l'entrée et de filtrer les signaux en est l'équivalent. Toutefois, lorsque vous utilisez la transformée de Fourier Discrète (DFT) (mis en œuvre avec un algorithme de transformée de Fourier Rapide pour la vitesse), vous avez réellement calculer un échantillonnage de la version du spectre. Cela a beaucoup de conséquences, mais l'une des plus pertinentes pour le filtrage est l'implication que le domaine temporel du signal est périodique.

Voici un exemple. Considérons un sinusoïdale du signal d'entrée x avec 1,5 cycles dans la période, et un simple filtre passe-bas h. Dans Matlab/Octave syntaxe:

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 glitch au début du bloc n'est pas ce que nous attendons tous. Mais si vous considérez fft(x), c'est logique. La transformée de Fourier Discrète suppose que le signal est périodique dans le bloc de transformation. Aussi loin que la DFT sait, nous avons demandé pour la transformation d'une période de ceci:

Cela conduit à la première considération importante lors de filtrage avec DFTs: vous êtes en fait la mise en œuvre de la convolution circulaire, pas de convolution linéaire. Donc le "glitch" dans le premier graphique n'est pas vraiment un problème si l'on considère les mathématiques. Alors la question devient: est-il un moyen de contourner la périodicité? La réponse est oui: l'utilisation overlap-save traitement. Essentiellement, vous devez calculer N-produits longs comme ci-dessus, mais seulement de 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 fait sens - nous nous attendre à un démarrage transitoire du filtre, le résultat s'installe dans l'état d'équilibre sinusoïdale de réponse. Notez que maintenant x peut être arbitrairement longue. La limitation est Nproc > 2*min(length(x),length(h)).

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

Il est difficile à voir, mais il y a un tas de non-zéro points à l'extrême bord gauche du graphique, puis un groupe de plus à l'extrême droite. En utilisant l'Octave intégré dans freqz fonction de regarder la réponse en fréquence nous voir (en faisant freqz(hraw)):

L'ampleur de la réponse a beaucoup de rides de la passe-haut de l'enveloppe vers le bas à zéro. Encore une fois, la périodicité inhérents à la DFT est à l'œuvre. Aussi loin que la DFT est concerné, hraw se répète encore et encore. Mais si vous prenez une période de hraw, comme freqz n', son spectre est assez différent de l'périodique de la version.

Donc, nous allons définir un nouveau signal: hrot = [hraw(513:end) ; hraw(1:512)]; Nous il suffit de tourner le raw DFT de sortie à faire en continu à l'intérieur du bloc. Maintenant, regardons la réponse en fréquence à l'aide de freqz(hrot):

Beaucoup mieux. L'souhaité enveloppe est là, sans ondulations. Bien sûr, la mise en œuvre n'est pas si simple maintenant, vous devez faire un complexe de multiplier par fft(hrot) plutôt que de simplement mise à l'échelle de chaque complexe bin, mais au moins vous aurez le droit de réponse.

Noter que pour la vitesse, vous feriez habituellement pré-calculer la TFD du collier de h, je l'ai laissé seul dans la boucle de comparer plus facilement avec l'original.

Êtes-vous d'atténuer la valeur de la DC fréquence de l'échantillon à zéro? Il semble que vous n'êtes pas d'en atténuer les effets dans votre exemple. Puisque vous êtes à la mise en œuvre d'un filtre passe-haut, vous devez définir la valeur DC à zéro en tant que bien.

Ceci pourrait expliquer la faible fréquence distorsion. Vous avez beaucoup de l'ondulation dans la bande de fréquence de réponse dans les basses fréquences, si le contrôleur de valeur est différente de zéro en raison de la grande transition.

Voici un exemple dans MATLAB/Octave pour démontrer ce qui se passe:

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

Avis que, dans mon code, je suis entrain de créer un exemple de la valeur DC non-zéro, suivi par un changement brusque de zéro, puis une rampe vers le haut. Je puis prendre la IFFT à transformer dans le domaine temporel. Puis-je effectuer un zéro-collier de fft (ce qui est fait automatiquement par MATLAB quand vous passez par une fft de taille plus grande que le signal d'entrée) sur le domaine temporel du signal. Le zéro-padding dans le domaine du temps, les résultats de l'interpolation dans le domaine des fréquences. Avec cela, nous pouvons voir comment le filtre répondra entre les échantillons du filtre.

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

Maintenant pour répondre à votre question sur la façon dont ce filtrage doit être fait. Il y a un certain nombre de façons, mais la plus simple à mettre en œuvre et de comprendre, c'est la conception de la fenêtre de la méthode. Le problème avec votre conception actuelle est que la largeur de la transition est énorme. La plupart du temps, vous voudrez aussi rapide de transitions que possible, avec aussi peu d'ondulation que 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'); 

Avis qu'il y a beaucoup d'oscillation provoquée par les changements brutaux.

La FFT ou de 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 -l'infini à l'infini, tandis que la DFT est appliquée sur un nombre fini d'échantillons. Cette en effet les résultats dans un carré de fenêtrage (troncature) dans le domaine temporel lors de l'utilisation de la DFT depuis nous ne nous occupons que d'un nombre fini d'échantillons. Malheureusement, la TFD d'un signal carré est une sinc fonction du type de (sin(x)/x).

Le problème avec le fait d'avoir des transitions nettes dans votre filtre (passe de 0 à 1 dans un échantillon), c'est que cela a une très longue réponse 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 du domaine temporel par une plus progressive de la fenêtre. Si l'on multiplie une fenêtre de hanning en ajoutant la ligne:

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

après la prise de la IFFT, on obtient cette réponse:

Donc je vous conseille d'essayer de mettre en œuvre la conception de la fenêtre de la méthode, car il est assez simple (il existe de meilleures façons, mais ils obtiennent plus compliqué). Puisque vous êtes à la mise en œuvre d'un égaliseur, je suppose que vous voulez être en mesure de modifier les atténuations à la volée, je vous suggère de calculer et de stocker le filtre dans le domaine fréquentiel à chaque fois qu'il y a un changement dans les paramètres, et alors vous pouvez simplement l'appliquer à chaque entrée de la mémoire tampon audio par la prise de la fft de la mémoire tampon d'entrée, en multipliant par votre domaine de fréquence du filtre échantillons, et en effectuant ensuite la ifft pour revenir à l'époque de domaine. Ce sera beaucoup plus efficace que toutes les branches que vous faites pour chaque échantillon.