Filterung - Multiplikation im Frequenzbereich

Ich versuche, ein einfaches Tiefpassfilter zu erstellen, aber ich habe das, was ich für ein überraschendes Ergebnis halte, wenn ich den Frequenzgang eines einfachen Butterworth-Filters betrachte.

Ich habe einen Großteil des folgenden Beispiels aus diesem anderen Beitrag kopiert . Ich habe am Ende des Skripts Code hinzugefügt, um die Eingangs- und Ausgangsspektren mit dem Frequenzgang des Filters zu vergleichen. Ich würde erwarten, dass das Ausgangsspektrum das Produkt des Eingangsspektrums und des Frequenzgangs : $\mathbf B$$\mathbf A$$\mathbf H$

$$\mathbf B = \mathbf H \mathbf A$$

Das Diagramm unten zeigt jedoch, dass der Filter tatsächlich einige Niederfrequenzkomponenten erhöht - sehen Sie, wie die rote Linie über der grünen unter .$4\textrm{ Hz}$

Könnte jemand bitte erklären, warum das so ist?

import numpy as np
from scipy.signal import butter, lfilter, freqz
import matplotlib.pyplot as plt
from scipy.fftpack import fft as fft
def butter_lowpass(cutoff, fs, order=5):
    nyq = 0.5 * fs
    normal_cutoff = cutoff / nyq
    b, a = butter(order, normal_cutoff, btype='low', analog=False)
    return b, a

def butter_lowpass_filter(data, cutoff, fs, order=5):
    b, a = butter_lowpass(cutoff, fs, order=order)
    y = lfilter(b, a, data)
    return y


# Filter requirements.
order = 6
fs = 30.0       # sample rate, Hz
cutoff = 3.667  # desired cutoff frequency of the filter, Hz

# Get the filter coefficients so we can check its frequency response.
b, a = butter_lowpass(cutoff, fs, order)

# Plot the frequency response.
w, h = freqz(b, a, worN=8000)
plt.subplot(2, 1, 1)
plt.plot(0.5*fs*w/np.pi, np.abs(h), 'b')
plt.plot(cutoff, 0.5*np.sqrt(2), 'ko')
plt.axvline(cutoff, color='k')
plt.xlim(0, 0.5*fs)
plt.title("Lowpass Filter Frequency Response")
plt.xlabel('Frequency [Hz]')
plt.grid()


# Demonstrate the use of the filter.
# First make some data to be filtered.
T = 5.0         # seconds
n = int(T * fs) # total number of samples
t = np.linspace(0, T, n, endpoint=False)
# "Noisy" data.  We want to recover the 1.2 Hz signal from this.
data = np.sin(1.2*2*np.pi*t) + 1.5*np.cos(9*2*np.pi*t) + 0.5*np.sin(12.0*2*np.pi*t)

# Filter the data, and plot both the original and filtered signals.
y = butter_lowpass_filter(data, cutoff, fs, order)

plt.subplot(2, 1, 2)
plt.plot(t, data, 'b-', label='data')
plt.plot(t, y, 'g-', linewidth=2, label='filtered data')
plt.xlabel('Time [sec]')
plt.grid()
plt.legend()

plt.subplots_adjust(hspace=0.35)
plt.show()

def calculateFFT(time,signal):
    N=len(signal)      
    df=1/((time[-1]-time[0]))
    frequencies=[i*df for i in range(int(N/2.0))]
    fftValues = [2.0/N*abs(i) for i in fft(signal,N)[0:N/2.0] ]
    return frequencies,fftValues

plt.subplot(2, 1, 1)
originalfreqs,originalFFT=calculateFFT(t,data)
plt.plot(originalfreqs,originalFFT,"g",label="original")

filteredfreqs,filteredFFT=calculateFFT(t,y)
plt.plot(filteredfreqs,filteredFFT,"r",label="filtered")
plt.legend()

Ich stimme dieser Frage gerne zu, denn so sollte eine gute Frage aussehen und auch, wie ich von den Schülern erwarten würde, dass sie ihr Verständnis für das, was sie lernen, überprüfen. Es ist immer gut, den Hintergrund von etwas zu verstehen, das Sie später verwenden möchten.

Das Problem, das Sie haben, ist das folgende: Im Prinzip haben Sie Recht, dass der Faltungssatz impliziert, dass die Faltung in einer Domäne eine Multiplikation in der anderen Domäne ist. Daher hast du

$$\mathcal{F}\{x(t)*h(t)\}=X(f)H(f)$$ wo $\mathcal{F}$ bezeichnet die Fourier-Transformation.

Also was ist $x(t)$ und $h(t)$ in Ihrem System? $x(t)$ist Ihr Eingangssignal (in blau), das die Summe von drei Sinus multipliziert mit einem rechteckigen Fenster ist. Es wird implizit mit einem Rechteck multipliziert, da Ihre Zeit nicht (in der Simulation nicht) reicht$-\infty$ zu $\infty$. Also das Spektrum von$X(f)$ist die Faltung einer sinc-Funktion (aus dem rechteckigen Fenster) mit drei Diracs (mit der Frequenz der Sinus). Dies ist eindeutig kein rein diskretes Spektrum.

Was ist $h(t)$? Es ist die Impulsantwort des Butterworth-Filters. Der Butterworth-Filter hat eine unendlich lange Impulsantwort, daher Ihr Faltungsprodukt$x(t)*h(t)$ist auch unendlich breit. Im Prinzip können Sie also keine endliche (dh diskrete) Fourier-Transformation (fft) anwenden und erwarten, dass sie dem kontinuierlichen Fall ähnlich ist.

Was Sie sehen, ist also vernünftig. Sie können versuchen, das Eingangssignal auf Null zu setzen, sodass Sie (den Hauptteil) der Impulsantwort in Ihrem Signal erhalten. Aber selbst dann werden Sie die Frequenzen sehen, da$X(f)$ ist eigentlich nicht diskret.

Zur Erweiterung auf @ Maximilian Matthé ‚s Antwort , können Sie die Visualisierung Spektralleckage Effekte (Faltung von sinc - Funktion in der Frequenzdomäne) durch Null-Auffüllen der Eingänge in calculateFFT. Mit der folgenden Funktion können Sie beispielsweise die Eingänge auf eine Länge auf Null setzen$k$ mal länger als die ursprüngliche Eingabe (wo in diesem Fall $k=4$):

def calculateFFT(time,signal):
    k=4
    N=k*len(signal)      
    df=1/(k*(time[-1]-time[0]))
    frequencies=[i*df for i in range(int(N/2.0))]
    fftValues = [k*2.0/N*abs(i) for i in fft(signal,N)[0:int(N/2.0)] ]
    return frequencies,fftValues

Wenn Sie das Ergebnis zeichnen, sollten Sie sehen, dass sich die Niederfrequenzkomponenten tatsächlich überlappen (was darauf hinweist, dass die Filterung das Signal dort nicht wirklich erhöht):

Woher kamen diese Wellen um den Spike? Es stellte sich heraus, dass sie immer da waren, aber bei der Berechnung der FFT erhielten Sie den Wert des Spektrums bei einem diskreten Satz von Frequenzwerten, und diese Welligkeiten passierten bei diesen Frequenzen genau Null. Wenn Sie eine etwas andere Signalfrequenz gewählt hätten, die kein genaues Vielfaches der FFT-Frequenzauflösung ist (in Ihrem Fall 0,2 Hz), z. B. 1,25 Hz anstelle von 1,2 Hz, hätte die Abtastung des Frequenzspektrums ganz anders ausgesehen (aufgrund zur Frequenzabtastung an verschiedenen Stellen in der Schwingung der Wellen):