DIY Fensterfunktionen

In der digitalen Signalverarbeitung werden Fensterfunktionen häufig verwendet, um das Signal zeitlich zu begrenzen, und ihre Namen sind jedem bekannt, der auf die eine oder andere Weise auf die diskrete Fourier-Transformation gestoßen ist: Hann, Hamming, Blackman, Harris und andere. Aber reichen sie aus, ist es möglich, sich etwas Neues auszudenken, und hat das einen Sinn?



In diesem Artikel betrachten wir die Ausgabe einer Fensterfunktion mit neuen Eigenschaften unter Verwendung von Wolfram Mathematica. Es wird auch davon ausgegangen, dass der Leser ein allgemeines Verständnis der digitalen Signalverarbeitung im Kontext des zur Diskussion stehenden Themas hat und zumindest mit dem Artikel aus Wikipedia vertraut ist .

Einführung

Historisch gesehen tauchten die ersten Fensterfunktionen bei Versuchen auf, ihre spektralen Eigenschaften durch Verringern der Amplitude der Nebenkeulen zu verbessern. Wenn das Signal mit der Fensterfunktion multipliziert wird, werden ihre Spektren gefaltet, was einige Einschränkungen für die Analyse mit sich bringt. Zu dieser Zeit, in den 50er Jahren des 20. Jahrhunderts, erlaubte die Rechenleistung nicht, symbolische Berechnungen einfach und natürlich zu manipulieren, und die Auswahl der optimalen Parameter war ziemlich schwierig. Dies ist einer der Gründe, warum die Funktionen mit unterschiedlichen Namen bezeichnet werden - sie, oder vielmehr ihre spektralen Eigenschaften, wurden von verschiedenen Forschern seit geraumer Zeit untersucht. Ein Nebeneffekt davon ist, dass der vorhandene Satz benannter Fensterfunktionen als eine Art "Standardsatz" wahrgenommen wird, außerhalb dessen nichts gefunden werden kann;Gleichzeitig enthalten die Namen dieser Fenster keine Informationen über die Eigenschaften und deuten auf eine separate Untersuchung und Speicherung hin.

Einstufung

Alle Fensterfunktionen werden erhalten, indem eine Funktion mit einer rechteckigen multipliziert wird, wodurch die Zeitbeschränkung garantiert wird. Daher können sie zunächst nach den von ihnen verwendeten Funktionen klassifiziert werden:

1. die Summe der Kosinusse. Die umfangreichste Klasse aufgrund der Tatsache, dass ihr Spektrum leicht berechenbar ist und die Summe der gewichteten Sinc-Funktionen ist. Dazu gehören Hann, Hamming, Blackman, Harris ...
2. stückweise kontinuierliches Polynom. Sie ergeben sich aus der Faltung einfacher Funktionen - zum Beispiel rechteckig und dreieckig. Gleichzeitig werden ihre Spektren multipliziert und ihr Auffinden ist auch nicht besonders schwierig. Dazu gehören Dirichlet, Triangular, Parzen, Welch ...
3. alle anderen - unter Verwendung von Exponentialen, Gaußschen, Sinc und anderen, ist die Wahl spezifischer Funktionen eher ideologischer Natur als spezifisch spektraler Eigenschaften.

Sie können auch durch Eigenschaften an den Rändern unterteilt werden:

1. Keine Lücken an den Kanten - Werte gleich Null. Die Pausen sind von Dirichlet, Hamming, Blackman-Harris. Nicht haben - Dreieckig, Hann, Nuttal
2. Fehlen von Diskontinuitäten des 1. und anderer Derivate

Zwei weitere Fensterfunktionen sind separat hervorzuheben:

1. Kaiser, mit dem Sie die Breite des Hauptblatts einstellen können;
2. Dolph-Chebyshev, dessen Nebenkeulen einer bestimmten Amplitude entsprechen.

Beide haben Diskontinuitäten an den Rändern von Werten und Ableitungen, und ihre Berechnungen sind mit einigen Schwierigkeiten verbunden - die Kaiser-Funktion wird durch eine spezielle Funktion (Bessel) berechnet, und die Dolph-Chebyshev-Funktion wird im Frequenzbereich bestimmt und durch die inverse diskrete Fourier-Transformation berechnet. Besonders schwierig ist es auch, ihre analytischen Antiderivate zu finden.

Reverse Engineering

Schauen wir uns die Definitionen der Blackman- und Nuttel-Funktionen an:

BlackmanWindow[x] // FunctionExpand

↓

$$

$$\left\{ \begin{array}{ll} \frac{1}{50} (25 \cos (2 \pi x)+4 \cos (4 \pi x)+21) & -\frac{1}{2}\leq x\leq \frac{1}{2} \\ 0 & \left| x\right| >\frac{1}{2} \\ \end{array} \right.$$

NuttallWindow[x] // FunctionExpand

↓

$$

$$\left\{ \begin{array}{ll} \frac{121849 \cos (2 \pi x)+36058 \cos (4 \pi x)+3151 \cos (6 \pi x)+88942}{250000} & -\frac{1}{2}\leq x\leq \frac{1}{2} \\ 0 & \left| x\right| >\frac{1}{2} \\ \end{array} \right.$$

Sie sind Summen von 3 und 4 Kosinuswellen mit geraden Frequenzen (beginnend bei Null) und einigen Koeffizienten. Woher kommen diese Koeffizienten? Es ist unwahrscheinlich, dass sie von Hand ausgewählt wurden, zumindest jeweils einzeln - schließlich gibt es Randbedingungen, die die Fenster unabhängig von ihren spektralen Eigenschaften erfüllen müssen - mindestens gleich eins in der Mitte.



Versuchen wir, die Blackman-Funktion selbst zu "erfinden". Dazu definieren wir eine Funktion mit noch unbekannten Koeffizienten

f = Function[x, a + b Cos[2 Pi x] + c Cos[4 Pi x]];

und ein Gleichungssystem zusammenstellen, das die Randbedingungen definiert - Gleichheit mit der Einheit in der Mitte, Null an den Rändern und Gleichheit mit Null der Ableitungen an den Rändern:

ss = Solve[{
   f[0] == 1,
   f[1/2] == 0,
   f'[1/2] == 0
   }, {a, b, c}]

↓

$$

$$\left\{\left\{b\to \frac{1}{2},c\to \frac{1}{2}-a\right\}\right\}$$

Es wurde die Lösung gefunden, aber nicht eine, sondern viele - abhängig vom Koeffizienten a , vor dem Wolfram uns höflich gewarnt hat. Aus der gefundenen Lösung definieren wir nun eine neue Funktion in Abhängigkeit vom unbekannten Koeffizienten: ↓ Es ist leicht zu erkennen, dass wir für a = 0,42 die Blackman-Funktion erhalten. Aber warum genau 0,42? Um diese Frage zu beantworten, müssen wir ihr Spektrum zeichnen. Die analytische Fourier-Transformation ist nicht die einfachste Aufgabe, aber Wolfram kann sie auch bewältigen und erspart uns so die Arbeit.



hx = Function[{x, a},

Piecewise[{{f[x] /. ss[[1]], -1/2 < x < 1/2}}, 0] // Evaluate]

hw = Function[{w, a}, 
  FourierTransform[hx[x, a], x, w] // #/Limit[#, w -> 0] & // 
    Simplify // Evaluate]

↓









In dieser Form ist der Graph der Funktion immer noch nicht sehr informativ, da es bequemer ist, das Spektrum auf einer logarithmischen Skala zu analysieren. Durch Hinzufügen der entsprechenden Umrechnung in Dezibel$$$20 \log _{10}(\left| x\right| )$erhalten wir den gleichen üblichen Graphen des Spektrums der Fensterfunktionen:









Beim Ändern des Parameters werden wir etwas Ähnliches beobachten:







Abhängig vom Parameter a ändert sich der Pegel der Nebenkeulen und bei a = 0,42 ist er mehr oder weniger minimal und gleichmäßig. Mit a = 0,409 können wir ein etwas besseres Ergebnis erzielen, wenn „besser“ als das minimal mögliche Niveau der Nebenkeulen verstanden wird.

Entwicklung

Offensichtlich war eine solche Mikrooptimierung die Mühe überhaupt nicht wert. Ist es möglich, die Eigenschaften dieses Fensters qualitativ zu verbessern, ohne die Anfangsdaten zu ändern?



Zuvor haben wir uns die Ausgabe von Fensterfunktionen angesehen, die sich zu eins summieren, sodass Sie das ursprüngliche Signal wiederherstellen können. Probieren wir diese Technik für unser Fenster aus und sehen, wie sie sich auf das Spektrum auswirkt.



Definieren wir eine Hilfsfunktion zum Erstellen überlappender Fenster unter Berücksichtigung der Grenzen der Argumentänderung im Bereich von -1/2 bis 1/2:







Suchen Sie das Antiderivativ, verschieben Sie es in die Mitte der Koordinaten und skalieren Sie es auf eins:





↓





Zeigen Sie es in der Grafik an:







Wie Sie sehen, hat Wolfram es auch hier alleine gemacht, und wir mussten keine stückweise kontinuierliche Definition eines Antiderivativs manuell angeben. Jetzt hängt die Ansicht unseres Fensters nicht nur von der Variablen a ab , sondern auch vom Grad der Überlappung - und wenn sie zunimmt, tendiert sie zur Form der Ableitung:







Der letzte Schliff besteht darin, eine analytische Funktion für das Spektrum zu finden, um den optimalen Wert des Parameters a zu bestimmen .



Wenn wir hier versuchen, die Transformation wie beim letzten Mal direkt zu berechnen, werden wir Wolfram in tiefe Gedanken treiben. Es gibt verschiedene Möglichkeiten, diesen Prozess zu beschleunigen:



- Verwenden Sie eine diskrete Transformation anstelle einer analytischen - die einfachste. Aber ein echter Mathematiker wird keine numerischen Methoden anwenden, bei denen eine rein analytische Lösung gefunden werden kann - und wir werden (noch) nicht; Darüber hinaus weist es offensichtliche Einschränkungen auf (in Bezug auf Auflösung und Spektrumbeschränkung).



- Übergeben Sie bestimmte Werte als Überlappungsparameter und versuchen Sie anhand der Ergebnisse, das Muster zu erkennen und zu verallgemeinern. Eine ziemlich funktionierende Option, die jedoch zusätzliche mentale und kreative Anstrengungen erfordert. Lassen wir es als letzten Ausweg.



- Berechnungen direkt im Frequenzbereich durchführen. Dies ist aufgrund der folgenden Eigenschaften der Fourier-Transformation möglich:



1) Sie ist linear, d.h. die Summe der Bilder ist gleich dem Bild der Summe;

2) Die Integration in den Zeitbereich entspricht der Division durch iw im Frequenzbereich.

3) Zeitversatz entspricht einer Modulation mit einer komplexen Sinuskurve.

(mehr auf Wikipedia oder hier ).



Wenn wir also das Spektrum hw einer beliebigen Fensterfunktion haben, können wir daraus das Spektrum hwo einer Funktion erhalten, die mit überlappendem t zu eins summiert wird :





↓





Nun können wir sehen, wie sich das Spektrum in der Dynamik ändert:







Hier wirkt sich eine Änderung des Überlappungsparameters bereits auf etwas andere Weise auf das Fensterspektrum aus :







Nachdem Sie ein wenig mit den Parametern gespielt haben, ist es leicht zu bemerken, dass "mehr nicht besser ist" und der optimale Grad der Überlappung für diese Funktion im Bereich von vier liegt. Insbesondere für t = 4 und a= 0,404 Wir erhalten einen Nebenkeulenpegel von nicht mehr als -80 dB. Dies ist ein sehr gutes Ergebnis - insbesondere wenn man bedenkt, dass die erhöhte Kosinusfunktion, die traditionell für summierte Fenster verwendet wird, einen Keulenpegel von etwa -30 dB ergibt. Wenn unsere neue Fensterfunktion explizit ausgeschrieben wird, sieht sie folgendermaßen aus:

$$

$$\left\{ \begin{array}{ll} \frac{404 \pi (1-2 x)+36 \sin \left(\frac{1}{3} (16 \pi x+\pi )\right)+375 \sin \left(\frac{1}{3} (\pi -8 \pi x)\right)}{606 \pi } & \frac{1}{4} < x \leq \frac{1}{2} \\ \frac{404 (2 \pi x+\pi )+36 \sin \left(\frac{2}{3} \pi (8 x+1)\right)+375 \sin \left(\frac{1}{3} (8 \pi x+\pi )\right)}{606 \pi } & -\frac{1}{2} < x \leq -\frac{1}{4} \\ \frac{\sqrt{3} \left(125 \cos \left(\frac{8 \pi x}{3}\right)+12 \cos \left(\frac{16 \pi x}{3}\right)\right)}{202 \pi }+\frac{1}{3} & -\frac{1}{4} < x \leq \frac{1}{4} \\ 0 & \left| x\right| >\frac{1}{2} \\ \end{array} \right.$$

Weitere Entwicklung

Was können Sie noch tun, um den Nebenkeulenpegel weiter zu senken? Sie können die Kosinusse nicht mit geraden, sondern mit ungeraden Frequenzen nehmen (hier ist aus Gründen der Kompaktheit die Lösung des Gleichungssystems direkt in die Definition der Funktion eingebettet):

$$

$$\left(a \cos (\pi x)+\left(\frac{5}{8}-\frac{a}{2}\right) \cos (3 \pi x)+\left(\frac{3}{8}-\frac{a}{2}\right) \cos (5 \pi x)\right) \Pi (x)$$

und nach Integration mit den Parametern a = 0,6628 und einem Überlappungspegel von 4,5 erhalten Sie eine Dämpfung von -90 dB:







Explizit:

$$

$$\small\left\{ \begin{array}{ll} \frac{327 \sin \left(\frac{1}{7} \pi (45 x+2)\right)+24855 \sin \left(\frac{1}{7} \pi (1-9 x)\right)+3670 \cos \left(\frac{1}{14} \pi (54 x+1)\right)+21512}{43024} & \frac{5}{18} < x \leq \frac{1}{2} \\ \frac{24855 \sin \left(\frac{1}{7} \pi (9 x+1)\right)+3670 \sin \left(\frac{3}{7} \pi (9 x+1)\right)+327 \sin \left(\frac{5}{7} \pi (9 x+1)\right)+21512}{43024} & -\frac{1}{2} < x \leq -\frac{5}{18} \\ \begin{array}{ll} \frac{-3670 \sin \left(\frac{3}{7} \pi (9 x-1)\right)+24855 \sin \left(\frac{1}{7} \pi (9 x+1)\right)+327 \sin \left(\frac{5}{7} \pi (9 x+1)\right)}{43024}+ \\ +\frac{3670 \sin \left(\frac{3}{7} \pi (9 x+1)\right)+327 \sin \left(\frac{1}{7} \pi (45 x+2)\right)+24855 \sin \left(\frac{1}{7} \pi (1-9 x)\right)}{43024} \end{array} & -\frac{5}{18} < x \leq \frac{5}{18} \\ 0 & \left| x\right| >\frac{1}{2} \\ \end{array} \right.$$

Sie können einen weiteren Kosinus hinzufügen und die Anzahl der Nullableitungen erhöhen:

$$

$$\left(a \cos (\pi x)+\frac{1}{80} (35-16 a) \cos (3 \pi x)+\frac{1}{80} (35-48 a) \cos (5 \pi x)+\frac{1}{40} (5-8 a) \cos (7 \pi x)\right) \Pi (x)$$

und nach Integration mit den Parametern a = 0,5862 und dem Überlappungspegel von 6,4 erhalten Sie eine Unterdrückung von -110 dB:







Explizit:

$$

$$\small\left\{ \begin{array}{ll} \frac{-3077550 \sin \left(\frac{1}{54} \pi (64 x-5)\right)-90069 \sin \left(\frac{5}{54} \pi (64 x-5)\right)-5820 \sin \left(\frac{7}{54} \pi (64 x-5)\right)-560455 \cos \left(\frac{1}{9} \pi (7-32 x)\right)+2601344}{5202688} & \frac{11}{32}< x \leq \frac{1}{2} \\ \frac{560455 \sin \left(\frac{1}{18} \pi (64 x+5)\right)+3077550 \sin \left(\frac{1}{54} \pi (64 x+5)\right)+90069 \sin \left(\frac{5}{54} \pi (64 x+5)\right)+5820 \sin \left(\frac{7}{54} \pi (64 x+5)\right)+2601344}{5202688} & -\frac{1}{2}< x \leq -\frac{11}{32} \\ \begin{array}{ll} \frac{-3077550 \sin \left(\frac{1}{54} \pi (64 x-5)\right)-90069 \sin \left(\frac{5}{54} \pi (64 x-5)\right)-5820 \sin \left(\frac{7}{54} \pi (64 x-5)\right)-560455 \cos \left(\frac{1}{9} \pi (7-32 x)\right)}{5202688}+ \\ +\frac{560455 \sin \left(\frac{1}{18} \pi (64 x+5)\right)+3077550 \sin \left(\frac{1}{54} \pi (64 x+5)\right)+90069 \sin \left(\frac{5}{54} \pi (64 x+5)\right)+5820 \sin \left(\frac{7}{54} \pi (64 x+5)\right)}{5202688} \end{array} &-\frac{11}{32}< x \leq \frac{11}{32} \\ 0 & \left| x\right| >\frac{1}{2} \\ \end{array} \right.$$

Eine noch signifikantere Verringerung des Niveaus der Nebenkeulen kann erreicht werden, indem die Fourier-Transformation quadriert wird und dadurch das Niveau der Nebenkeulen gleichzeitig um den Faktor 2 verringert wird. Dies ermöglicht es Ihnen, die Zunahme der Anzahl von Parametern für ihre manuelle Auswahl zu beseitigen, erhöht jedoch die Komplexität bei der Berechnung der Faltung im Zeitbereich.









und hier ist es bereits möglich, eine Unterdrückung über 160 dB zu erreichen:







Wir werden die Formel wegen ihrer beeindruckenden Größe nicht explizit angeben.

Hypergeometrische Fensterfunktionen

Um die erforderliche Anzahl von Nullen an den Grenzen unserer Funktion sicherzustellen, haben wir die Lösung eines Gleichungssystems mit anschließender Integration durchsucht. Dies ist nicht sehr praktisch - Sie müssen die Anzahl der Gleichungen jedes Mal ändern. Vielleicht gibt es eine einfachere und schönere Lösung? Es gibt! Und die hypergeometrische Funktion hilft uns dabei.





Unsere Lösung sieht folgendermaßen aus:

$$

$$\frac{2 z \Gamma \left(\frac{n+3}{2}\right) \, _2F_1\left(\frac{1}{2},-\frac{n}{2};\frac{3}{2};z^2\right)}{\sqrt{\pi } \Gamma \left(\frac{n}{2}+1\right)}+a \left(z+b z^3\right) \left(1-z^2\right)^{n/2}$$

Wo $$$z=\sin \left(\frac{\pi x}{2}\right)$



Es besteht aus der Summe zweier Teile: dem Antiderivativ der Funktion $$$\left(1-z^2\right)^{\frac{n}{2}}$, in der Amplitude auf den Punkt (1,1) normiert und eine Funktion mit freien Parametern, multipliziert mit derselben $$$\left(1-z^2\right)^{\frac{n}{2}}$Um die erforderliche Anzahl von Nullableitungen beizubehalten: Die







Genauigkeit der Anpassung und ihr Einfluss auf die Verteilung der Nebenkeulen hängt davon ab, welche Funktion wir zur Korrektur der Form wählen. Es gibt mögliche Optionen, die separate Studien erfordern. In diesem Fall die Funktion$$$a \left(z+b z^3\right)$bietet eine völlig akzeptable Option - und aufgrund von zwei Parametern können Sie genauere Einstellungen vornehmen.



Die Formel selbst ist insofern interessant (und praktisch), als sie für gerade n zu einem Polynom und für ungerade n zu der Summe eines Polynoms mit Arkussinus und Quadratwurzel vereinfacht wird - was die endgültige Formel kompakter und einfacher zu berechnen macht.



Es ist bereits einfacher (und schneller), die Parameter hier durch die diskrete Fourier-Transformation auszuwählen. Wir benötigen außerdem eine zusätzliche Definition der Überlappungsfunktion, um mit drei Parametern arbeiten zu können:



















Als Beispiel wird unsere Funktion nach dem Ersetzen der Parameter aus dem Bild vereinfacht

$$

$$\frac{-288 z^{11}+2315 z^9-7380 z^7+12330 z^5-11940 z^3+8163 z}{3200}$$

Inverses Kaiserfenster

Wenn sich die Aufgabe, die Fenster zu einem zusammenzufassen, nicht lohnt, ist das Kaiser-Fenster das optimalste Fenster, mit dem Sie die Breite der Hauptkeule reibungslos anpassen können. Es hat jedoch einen Nachteil - da es in Form der Bessel-Funktion ausgedrückt wird , ist es etwas schwierig, es außerhalb der mathematischen Pakete zu lesen. Sie können die Bessel-Funktion natürlich separat implementieren - oder Sie können nach einer Annäherung durch Elementarfunktionen suchen. Und plötzlich stellte sich heraus, dass man mit der Funktion des zeitlich begrenzten Spektrums des Kaiser-Fensters ein noch etwas besseres Ergebnis erzielen kann -

$$

$$\frac{ \sinh \left(k \sqrt{1-4 x^2}\right)}{\sinh (k) \sqrt{1-4 x^2} } \Pi (x)$$

Die Bezahlung für eine solche listige Entscheidung war die Unmöglichkeit, ihr Spektrum rein analytisch auszudrücken - aber das ist nicht so wichtig, weil wir in der Praxis auf jeden Fall Daten in diskreter Form mit der diskreten Fourier-Transformation betreiben und analysieren. Unten in der Grafik können Sie ihre Spektren vergleichen, wobei Rot das ursprüngliche Kaiser-Fenster und Grün die Annäherung ist:









Bei gleicher Breite des Hauptblatts ist die Höhe der Nebenlappen etwas geringer. Zur Vereinfachung der Verwendung können Sie eine Tabelle mit ungefähren Werten des Parameters k erstellen , um das erforderliche Maß an Nebenkeulenunterdrückung zu erhalten:

$$

$$\begin{array}{cc} \text{-60 } & \text{k=8.8} \\ \text{-90 } & \text{k=11.36} \\ \text{-120 } & \text{k=15.18} \\ \text{-150 } & \text{k=18.88} \\ \end{array}$$

Das Finden einer Annäherung an das Antiderivativ des inversen Kaiserfensters bleibt ein separates interessantes Problem. Bisher war es nur möglich, es durch eine Potenzreihe auszudrücken, die ziemlich langsam konvergiert:

$$

$$\small \int \frac{ \sinh \left(k \sqrt{1-4 x^2}\right)}{\sinh (k) \sqrt{1-4 x^2} } \, dx= \sum _{n=1}^{\infty } -\frac{2^{\frac{1}{2}-n} k^{n-1} x^{2 n-1} \left(\sqrt{-k} (-1)^n K_{n-\frac{1}{2}}(-k)+\sqrt{k} K_{n-\frac{1}{2}}(k)\right)}{\sqrt{\pi } (2 n-1) \sinh (k) \Gamma (n)}$$

Vielleicht können Experten in Mathematik in den Kommentaren eine bessere Lösung vorschlagen.

Fazit

Trotz der Tatsache, dass die digitale Signalverarbeitung eine vollständig ausgebildete Disziplin ist, gibt es noch Raum für Entwicklung und etwas Neues darin - zumindest in der vorhandenen speziellen und wissenschaftlichen Literatur werden die Fragen der Konstruktion von Fensterfunktionen, die zu einer Einheit zusammengefasst werden, nicht berücksichtigt. und die Fähigkeiten moderner Computeralgebrasysteme ermöglichen es, das gesammelte Wissen auf einem qualitativ neuen Niveau zu überdenken.



Das Originaldokument von Wolfram Mathematica mit allen Berechnungen ist auf GitHub verfügbar .