Fouriertransformation über Teilintervall

Gregor

Hi.

Ich habe hier folgendes Problem:

Es ist eine periodische Funktion mit einer langen Periode über ein sehr großes Intervall gegeben. Die Funktion ist frequenzbeschränkt und sie liegt in Form von Koeffizienten einer Fouriertransformierten perfekt abgetastet vor.

Jetzt möchte ich etwas rechnerisch aufwändiges mit der Funktion machen, für das ich aber nur das Verhalten der Funktion in einem kleinen Teilbereich des Intervalls kennen muss. Ich muss es allerdings auch zwischen den Abtastpunkten kennen. Um den Rechenaufwand zu verringern möchte ich nun die Funktion in den Realraum transformieren, dort einen Bereich um das relevante Teilintervall herum auswählen und die Funktion über dieses Teilintevall wieder in den reziproken Raum zurücktransformieren. Ich hätte dann eine Funktion im reziproken Raum gegeben, die eine kürzere Periode hat und somit durch weniger Fourierkoeffizienten darstellbar ist.

Soweit nur die Funktionswerte an den Abtastpunkten relevant wären, wäre das soweit kein Problem. Das Problem ist jetzt, dass an den Grenzen des Teilintervalls die neue Funktion mit der kürzeren Periode Unstetigkeitsstellen hat. Dadurch ist die Funktion nicht mehr frequenzbeschränkt und man sieht zwischen den Abtaststellen Gibbs-Oszillationen. Die möchte ich aber vermeiden, weil ich das Verhalten der Funktion ja auch zwischen den Abtastpunkten brauche.

Eine Idee zur Vermeidung der Oszillationen wäre jetzt, das Teilintervall etwas größer als den für mich relevanten Bereich zu wählen und in dem Randbereich, der für mich nicht relevant ist, die Funktion durch einen künstlichen Funktionsverlauf zu ersetzen, der Stetigkeit (auch in den Ableitungen) an den Intervallgrenzen garantiert und die Maximalfrequenz der Funktion nicht erhöht.

Die Frage ist jetzt, wie man das besonders schlau macht. Ich gehe davon aus, dass so eine Problemstellung häufiger auftritt und, dass es deshalb Standardlösungen dafür gibt. Weiß jemand, wie man normalerweise mit so etwas umgeht? Oder gibt es einen prinzipiell anderen Ansatz, um mit so einem Problem umzugehen.

Vielen Dank schon mal für alle konstruktiven Ideen im Voraus!

nachtfeuer

Ich gehe mal davon aus, dass dir dieser Hintergrund sehr vertraut ist
-> https://de.wikipedia.org/wiki/Nyquist-Shannon-Abtasttheorem

Früher bei Synthesizern und Samplern war das in etwa so: man hatte nicht viel Speicher und sehr begrenzte Bitbreiten.
Bei Synthesizern nahm man sich, um Speicherplatz zu sparen, sehr kurze Samples (z.B. für Standardwellenformen) und Psychoakustik zu Hilfe.
So bekam man z.B. mit einem guten (und kurzen) Anblassample + synthetisierte Rechteckwellenform dazu einen gut klingenden Flötensound.
Manchmal waren die per MidiSampleDump erlaubten Längen (einiger Geräte) aber viel zu kurz, und es kam nicht unbedingt der gewünschte Sound (komplexere Wellenform) heraus (oft eher eine Art krankes Weckerrattern oder ähnlich).

Bei Samplern war es immer auch eine Kunst, einigermaßen brauchbare Loopstellen zu finden.

Die Ansätze da waren vor allem Filter, Hüllkurven/LFOs, Up- und Downsampling, Kombinationen (oder eben Stacks oder Wellenmix) und viel damit herumspielen und eben längere Samplezeiten und Standardwellenformen.
-----------------------------------------------

Was jetzt auch noch nicht klar ist, hast du eher nur ein mathematisches und/oder auch noch ein Compilerproblem/bzw. Digitalisierproblem dazu.

(Zusatzfrage) Was würde passieren, wenn du ein anderes Zahlensystem zugrunde legst (z.B. Hexadezimal, Festkomma, Fließkomma etc.)?

Gregor

nachtfeuer schrieb:

Was jetzt auch noch nicht klar ist, hast du eher nur ein mathematisches und/oder auch noch ein Compilerproblem/bzw. Digitalisierproblem dazu.

(Zusatzfrage) Was würde passieren, wenn du ein anderes Zahlensystem zugrunde legst (z.B. Hexadezimal, Festkomma, Fließkomma etc.)?

Das ist eine rein mathematische Fragestellung. Bei der Funktion geht es auch nicht um ein Tonsignal, sondern um eine Elektronendichte in einem Kristall. Ich brauche die Funktion zwischen den Abtaststellen auf mehrere (vielleicht 6 oder so) Stellen genau.

Ich glaube, ein Ansatz über eine andere Darstellung der Werte der Funktion wird mich nicht weiterbringen. Das löst das Problem auf der mathematischen Ebene nicht.

Jodocus

Hi!

Ich bin jetzt auch nicht direkt ein Spezialist für derart Probleme, aber bietet sich hier nicht evtl. ein Mollifier an, sodass du nach Faltung mit der Indikatorfunktion über dein Intervall eine glatte Bumpfunktion erhälst? Wenn du die mit deiner Funktion multiplizierst, hast du eine glatte Funktion, die in gewünschten Intervall exakt die gleichen Werte hat und innerhalb der beiden Randintervalle glatt auf 0 abfällt, sich also stetig/diffbar periodisch fortsetzen lässt. Inwieweit du dadurch dann mehr Fourierkomponenten bekommst bzw. welcher Mollifier am besten ist (im Sinne von, bringt möglichst wenig Fourierkomponenten ins Endresultat), kann ich dir jetzt pauschal auch nicht sagen...

Gregor

Jodocus schrieb:

Hi!

Ich bin jetzt auch nicht direkt ein Spezialist für derart Probleme, aber bietet sich hier nicht evtl. ein Mollifier an, sodass du nach Faltung mit der Indikatorfunktion über dein Intervall eine glatte Bumpfunktion erhälst? Wenn du die mit deiner Funktion multiplizierst, hast du eine glatte Funktion, die in gewünschten Intervall exakt die gleichen Werte hat und innerhalb der beiden Randintervalle glatt auf 0 abfällt, sich also stetig/diffbar periodisch fortsetzen lässt. Inwieweit du dadurch dann mehr Fourierkomponenten bekommst bzw. welcher Mollifier am besten ist (im Sinne von, bringt möglichst wenig Fourierkomponenten ins Endresultat), kann ich dir jetzt pauschal auch nicht sagen...

Ja, an so etwas hatte ich gedacht. Den Begriff Mollifier kannte ich noch nicht, obwohl ich ähnliches schon häufiger zur Approximation eines Delta-Impulses genutzt habe.

Jetzt ist nur noch die Frage, wie die perfekte Form des Mollifiers ist.

Was ist denn da genau der Unterschied zu https://de.wikipedia.org/wiki/Fensterfunktion ?

krümelkacker

Gregor schrieb:

Es ist eine periodische Funktion mit einer langen Periode über ein sehr großes Intervall gegeben. Die Funktion ist frequenzbeschränkt und sie liegt in Form von Koeffizienten einer Fouriertransformierten perfekt abgetastet vor.

Jetzt möchte ich etwas rechnerisch aufwändiges mit der Funktion machen, für das ich aber nur das Verhalten der Funktion in einem kleinen Teilbereich des Intervalls kennen muss. Ich muss es allerdings auch zwischen den Abtastpunkten kennen.

Wenn ich das richtig verstanden habe, willst Du auf einem Teilausschnitt der periodischen Funktion, die Du als Koeffizienten der Fourier-Transformation vorliegen hast, "arbeiten" und benötigst dabei "Zwischenwerte".

Ich würde das in zwei Schritten machen:

Eine Interpolationsstufe in der Frequenzdomäne (Fourier-Koeffizienten-Darstellung) per Zero-Padding
Eine zweite Interpolationsstufe in der Zeit/Ortsdomäne auf einem kleinen Teilausschnitt per Zero-Stuffing + Lowpass Filtering

Das kann dann so aussehen:

signal = randn(1,16);  % kritisch abgetastetes, zufälliges, periodisches Signal
coeffs = fft(signal);  % Frequenz-Darstellung

oversampling = 4;      % Überabtastung
n = length(coeffs);
padded = [coeffs(1:n/2) coeffs(n/2+1)*0.5 zeros(1,n*(oversampling-1)-1) coeffs(n/2+1)*0.5 coeffs(n/2+2:end)] * oversampling;
signal2 = real(ifft(padded)); % real() nur benutzen, wenn wirklich reellwertig!

% Plotten
hold off;
plot((0:n*oversampling)./(n*oversampling),[signal2 signal2(1)],'g');
hold on;
stem((0:n-1)./n,signal,'r');

% Zweite Interpolationsstufe um Faktor 8
oversampling2 = 8;
b = firls(5*oversampling2,[0 1/oversampling 2-1/oversampling oversampling2]./oversampling2,[1 1 0 0]);
signal3 = upfirdn(signal2(20:40)*oversampling2,b,oversampling2,1);

plot((132+(0:length(signal3)-1))./(n*oversampling*oversampling2),signal3,'k');

signal2 hat jetzt 4mal soviele Stützstellen wie signal und signal3 hat jetzt 8mal soviele Stützstellen wie ein Ausschnitt aus signal2 . b ist die Impulsantwort eines Tiefpassfilters, der für die zweite Interpolationsstufe benutzt wird. Dadurch, dass signal2 überabgetastet ist, kann die Impulsantwort b relativ kurz sein und dementsprechend braucht man auch nicht viele Werte außerhalb des Bereichs, der dich interessiert.

Hier der Plot des Skripts.

Gregor

Daniel@work schrieb:

Was ist denn da genau der Unterschied zu https://de.wikipedia.org/wiki/Fensterfunktion ?

Ein berechtigter Einwand. Im Prinzip läuft das bisher auf so etwas hinaus, wobei dann die Frage ist, welche Fensterfunktion die Beste ist, um sie mit der Indikatorfunktion zu falten. Aber auf der von Dir verlinkten Seite sind ja jede Menge Frequenzspektren zu sehen, so dass das eine gute Entscheidungshilfe wäre.

Eigentlich brauche ich aber mehr. Wir wissen ja, dass Funktionen maximal entweder im Realraum oder im Frequenzraum einen kompakten Träger haben können. Die Funktion, die ich konstruieren möchte, ist aber im Realraum periodisch und soll im Frequenzraum einen kompakten Träger haben. Ich glaube, ich kann mein Ziel durch die Faltung/Multiplikation mit einer Fensterfunktion nicht erreichen. Es sei denn, es gibt irgendwo das Konzept von periodischen Fensterfunktionen, die im Realraum kompakt sind. Habe da aber auf die Schnelle nichts zu gefunden.

@krümelkacker: Ja, die Idee, dass ich Zeropadding betreiben muss, ist mir auch schon gekommen. Letztendlich ist es ja so, dass wenn ich die Funktion mit irgendeiner anderen Funktion multipliziere, um sie an den Intervallrändern auf Null zu drücken, ich dann auf jeden Fall höhere Frequenzen hereinkriege. So eine Multiplikation im Realraum ist schließlich eine Faltung im Frequenzraum.

Ich denke momentan über eine völlig andere Herangehensweise nach:

Da ich von einer Funktion im Frequenzraum als Startpunkt ausgehe, kann ich von der Funktion im Realraum problemlos die exakten 1-ten und 2-ten Ableitungen an jedem Abtastpunkt bestimmen. Jetzt formuliere ich die Bandbeschränkung mal so, dass ich die "kinetische Energie" in den Frequenzen oberhalb der Grenzfrequenz auf Null drücken möchte. Die "kinetische Energie" kann ich sowohl im reziproken Raum als auch im Realraum berechnen. Im Realraum komme ich da über die 2-ten Ableitungen ran, im reziproken Raum kriegt man das direkt über die Funktionswerte zu jeder Frequenz.

Mir ist noch nicht klar, ob ich mein Ziel so umformulieren kann oder ob ich eigentlich mehr Bedingungen brauche.

Wenn ich das machen kann, dann könnte ich zumindest ein lineares Gleichungssystem mit den Funktionswerten und 1-ten und 2-ten Ableitungen aufstellen, die meine Linearkombination aus trigonometrischen Funktionen erfüllen muss. Ich könnte von Anfang an genau kontrollieren, wie viele zusätzliche Abtastpunkte ich in meinem Intervallrand benötige, weil ich genau wüsste, wie viele Funktionen ich brauche, um mein lineares Gleichungssystem lösen zu können.

Ich weiß, dass mal ein ähnlicher Ansatz in einem anderen Zusammenhang (Paywall) gemacht wurde. Hier wurde eine zu konstruierende Funktion durch eine Basis jenseits von trigonometrischen Funktionen ausgedrückt und ein Optimierungsproblem gelöst, um bei der Entwicklung in diese Basis die kinetische Energie jenseits einer Grenzfrequenz zu minimieren.

Sind Euch ähnliche Ansätze in solchen Zusammenhängen bekannt?

C14

Was optimal ist, hängt bestimmt davon ab, was du eigentlich am Ende des Tages berechnen willst.
Kannst du noch etwas mehr zum Problem sagen?
Das Problem ist eindimensional, richtig?
Wie groß ist das Teilinterval größenordnungsmäßig verglichen mit dem Gesamtinterval und der Abtastrate?
Welche Berechnung stellst du mit der Funktion in dem Teilintervall an?
(an welchen Stellen brauchst du dafür den Funktionswert?
Brauchst du den überhaupt oder vielleicht eher Skalarprodukte mit anderen Funktionen, die du evtl. direkter approximieren kannst?)

Gregor

C14 schrieb:

Was optimal ist, hängt bestimmt davon ab, was du eigentlich am Ende des Tages berechnen willst.
Kannst du noch etwas mehr zum Problem sagen?
Das Problem ist eindimensional, richtig?
Wie groß ist das Teilinterval größenordnungsmäßig verglichen mit dem Gesamtinterval und der Abtastrate?
Welche Berechnung stellst du mit der Funktion in dem Teilintervall an?
(an welchen Stellen brauchst du dafür den Funktionswert?
Brauchst du den überhaupt oder vielleicht eher Skalarprodukte mit anderen Funktionen, die du evtl. direkter approximieren kannst?)

Ich habe eigentlich versucht, das Gesamtproblem auf das Wesentliche zu reduzieren, aber hier noch etwas mehr dazu:

Das Problem ist eigentlich dreidimensional. Ich habe eine periodische Funktion und muss deren Werte und Ableitungen auf Kugeloberflächen kennen. Die (nicht-überlappenden) Kugeln sind hierbei im Realraum innerhalb einer repräsentierenden Periode der Funktion gegeben. Um die Werte und Ableitungen der Funktion auf der Kugeloberfläche zu bekommen, wird momentan eine Rayleigh-Entwicklung von ebenen Wellen nach Kugelflächenfunktionen durchgeführt.

Der Punkt ist jetzt die Zeitkomplexität bei der Skalierung des Problems. Die Skalierung des Problems ist im Wesentlichen eine Vergrößerung des periodisch fortgesetzten Bereichs. Hierbei wächst auch die Anzahl der Kugeln N_K innerhalb dieses Bereichs. Es gibt innerhalb des periodisch fortgesetzten Bereichs etwa 1000*N_K Abtastpunkte, bzw. die Anzahl der ebenen Wellen, mit denen ich die Funktion im Frequenzraum beschreibe, ist auch 1000*N_K.

Wenn ich nun für alle Kugeln diese Rayleigh-Entwicklung durchführe, dann skaliert das also wie 1000*N_K^2. Dieses quadratische Skalierungsverhalten ist eigentlich nicht wirklich zu rechtfertigen, weil ich durch die Vergrößerung des periodischen Bereichs lokal keine zusätzlichen Informationen über die Funktion erhalte und an jeder Kugel nur das lokale Verhalten der Funktion benötige. Das Ganze sollte eigentlich mehr oder weniger linear mit der Problemgröße skalieren.

Die Rayleigh-Entwicklung beruht darauf, dass man die Funktion als Linearkombination ebener Wellen gegeben hat. Der Ausgangspunkt für diese Entwicklung ist also eine Darstellung im Frequenzraum. Die Idee ist jetzt, dass ich die ursprüngliche Frequenzraumdarstellung mittels einer schnellen Fouriertransformation in O(N_K*log(N_K)) in den Realraum transformiere, dort um jede Kugeln herum in O(N_K) durch lokale Hilfsfunktionen ersetze und dann für jede Kugel die Hilfsfunktion in O(N_K) zurück in den Frequenzraum transformiere. Die Rayleigh-Entwicklungen würden dann zusammen in O(N_K) ablaufen, da der lokale Bereich für jede Kugel nicht mit der Problemgröße skaliert. Durch so einen Ansatz würde ich insgesamt also von einem quadratischen Skalierungsverhalten hin zu einem O(N_K*log(N_K)) Verhalten kommen, wobei es einen additiven linearen Term gibt, der die Zeitkomplexität für die relevanten Problemgrößen (100 < N_K < 10000++) dominieren sollte. Der O(N_K*log(N_K)) Beitrag zur Skalierung kommt nur dadurch, dass die Funktion ursprünglich in einer Frequenzraumdarstellung gegeben ist, davon wird man nicht wegkommen.

Ich brauche am Schluss keine Skalarprodukte, sondern wirklich das Verhalten der Funktion und ihrer Ableitungen auf der Kugeloberfläche. Punktweise und auf mehrere Stellen genau.