Grenzwerte sinc(x) bei float, bzw. double

Bashar

Ich würde erwarten, dass ein praktisch relevanter Faltungskern wesentlich schmaler ist als der, den du bekommst, wenn du nur die exakten Nullen an den Rändern weglässt.

(BTW grau ist alle Theorie:
sinc(1e300) = -2.8726e-301 != 0, damit hat sich die ursprüngliche Frage wohl erledigt)

Gibt es da womöglich eine praktische Obergrenze, unter der sinc(x) in der Berechnung keinen wirklichen Unterschied mehr macht? Wie gesagt, ab einer bestimmten Größe sind die Ergebnisse da eh mehr oder weniger zufällig, d.h. du kriegst dann etwas, was vermutlich auf leises Weißrauschen hinausläuft. Wenn x erstmal 10¹⁶ erreicht hat (jetzt mit ieee754-doubles gerechnet), springst du da in Ganzzahlschritten; spätestens dann verliert sinc(x) jede Aussagekraft.

Ansonsten würde ich mich da eher an dem bereits bestehenden Ergebnis orientieren. Du kannst eine Obergrenze für die Größenordnung von sinc(x) aus x abschätzen, und bei einer Summe

out[i] += sinc(x);

macht sinc(x) (wieder mit ieee754-doubles) dann keinen Unterschied mehr, wenn out[i] größer als etwa 10¹⁶ * sinc(x) ist. Für float und long double kannst du die Mantissenlänge aus std::numeric_limits ablesen; Bei ieee754-singles läuft das (glaube ich) auf etwa 8 Stellen hinaus, bei long doubles ist es unterschiedlich. Wenn du vorher abschätzen kannst, wie groß out[i] am Ende werden wird (vielleicht stellt sich das ja zu Anfang gleich heraus), kannst du auf die Art eine praktische Obergrenze für x festlegen, ab der sinc(x) < 1 / x in der Fließkommarechnung verschwindet. Das ginge auch on the fly, aber wenn du den Aufwand pro Schleifendurchlauf drücken willst, ist das womöglich kein gangbarer Weg.

Es ist aber, wie gesagt, davon auszugehen, dass du da am falschen Grenzwert forschst.

PhilippHToner

seldon schrieb:

dass das bei x > DBL_MAX / M_PI der Fall sein dürfte.

Da wir hier q/(x*M_PI) haben denke ich auch, dass, wenn der Nenner DBL_MAX ist, das Ergebnis auch 0 sein müssen, also für x=DBL_MAX/M_PI, aber ich hoffe, dass es eine frühere Grenze gibt.

krümelkacker schrieb:

Aber warum interessiert Dich das?

http://www.dspguide.com/ch16.htm

Die Möglichkeit den sinc zu fenstern kenn ich auch schon, allerdings würde ich gerne auch die fehlerfreieste möglich darstellbare Möglichkeit implementieren.

Bashar schrieb:

Ich würde erwarten, dass ein praktisch relevanter Faltungskern wesentlich schmaler ist als der, den du bekommst, wenn du nur die exakten Nullen an den Rändern weglässt.

(BTW grau ist alle Theorie:
sinc(1e300) = -2.8726e-301 != 0, damit hat sich die ursprüngliche Frage wohl erledigt)

Ja also ich brauche noch einige Samples mehr, dass der sinc 0 werden würde, danke daran habe ich gar nicht gedacht, einfach mal so auszuprobieren.

Jetzt aber mal eine andere Frage: Wann wird der sin zu ungenau? Ich hatte eigentlich gedacht, dass mir Google anwortet, aber dem ist nicht so.
Sollte ich im Sinus einen fmod einbauen?

sinc(x) = sin([e]pi[/e]*fmod(x, 2.f)) / ([e]pi[/e]*x)

fmod wird dir nicht weiterhelfen, weil sie in die selben Probleme läuft wie sin. Ein double x hat halt nur 53 signifikante binäre Stellen (was in der Praxis etwa 16 signifikante Dezimalstellen bedeutet). Wenn x > 10¹³ ist, hast du noch etwa drei Nachkommastellen. Wenn x > 10¹⁴, dann sind es noch zwei, und so weiter. Wie viele brauchst du denn, damit sinc(x) noch Sinn ergibt?

PhilippHToner

seldon schrieb:

Ansonsten würde ich mich da eher an dem bereits bestehenden Ergebnis orientieren. Du kannst eine Obergrenze für die Größenordnung von sinc(x) aus x abschätzen, und bei einer Summe
out[i] += sinc(x);
macht sinc(x) (wieder mit ieee754-doubles) dann keinen Unterschied mehr, wenn out[i] größer als etwa 10¹⁶ * sinc(x) ist. [...]

Es ist aber, wie gesagt, davon auszugehen, dass du da am falschen Grenzwert forschst.

Es ist eine Faltung. Der Sinc wird mit den Samples gewichtet. Dein Code würde das Sinc Integral berechnen, welches immer gleich ist. Hier mal mein jetziger funktionierender Code:

template<typename Sampling, typename Precision>
	void ResamplingSystem::perform(Sampling &newSamples) const
	{
#if _DEBUG
			std::cout << "Perform sampling..." << std::endl;
#endif
		const Precision samplesPerNewSample = (Precision)m_samples.size() / (Precision)newSamples.size();

		for(size_t index = 0; index < newSamples.size(); ++index)
		{
			Precision sample = (Precision)0;

			size_t oldIndex = m_samples.size();
			while(oldIndex--)
				sample += m_samples[oldIndex] * sinc<Precision>( samplesPerNewSample * (Precision)index - (Precision)oldIndex);

			newSamples[index] = sample;

#if _DEBUG
			std::cout << index << ", ";
#endif
		}
#if _DEBUG
		std::cout << std::endl << "Finished!" << std::endl;
#endif
	}

Ich denke ich werde sinc fenstern und die Fensterbreite variabel machen. Jenseits des Fenster ist dann alles 0, somit muss nicht vollständig gefaltet werden.
Mir fällt sowieso gerade ein, dass ich ohnehin den sinc fenstern muss, nachdem die Grenzen nicht 0 sind. Somit müsste ich ohnehin bei meinem Beispiel den sinc über die komplette Breite fenstern.

Danke nochmal für die Hilfen!

PhilippHToner

seldon schrieb:

fmod wird dir nicht weiterhelfen, weil sie in die selben Probleme läuft wie sin. Ein double x hat halt nur 53 signifikante binäre Stellen (was in der Praxis etwa 16 signifikante Dezimalstellen bedeutet). Wenn x > 10¹³ ist, hast du noch etwa drei Nachkommastellen. Wenn x > 10¹⁴, dann sind es noch zwei, und so weiter. Wie viele brauchst du denn, damit sinc(x) noch Sinn ergibt?

Ich hatte halt gedacht, wenn der sin(a) zuerst zu fehleranfällig wird, rechne ich halt mit sin( fmod(a, 2*M_PI)). Somit ist die math. Reihe doch genauer?

Das bringt dir doch nicht magisch mehr signifikante Stellen. Du hast beispielsweise (wir tun mal so, als wären Fließkommazahlen in Radix 10 gespeichert)

x           =  1234567890123.456
fmod(x, 1.) = 0.4560000000000000

Das ist natürlich in der Praxis binär und sieht nur in Binärdarstellung so aus, aber alles, was hier nach der dritten Nachkommastelle kommt, ist auch nach fmod reiner Zufall.

Wenn du x vorher in der Berechnung klein halten kannst, dann mag das was anderes sein, aber das kann ich dir ohne Kenntnis der Materie natürlich nicht sagen.

SeppJ

PhilippHToner schrieb:

Ich hatte halt gedacht, wenn der sin(a) zuerst zu fehleranfällig wird, rechne ich halt mit sin( fmod(a, 2*M_PI)). Somit ist die math. Reihe doch genauer?

Das Problem ist aber (auch), wie genau der Ursprungswert ist. Die Fließkommazahlen liegen bei großen Zahlen immer weiter auseinander. Wenn das a nur noch auf +-10 genau ist, dann nützt dir auch ein fmod nichts.

Tachyon

Erstmal eine Frage: Ist das Umtastverhätnis konstant?

Außerdem: Was Du vorhast, mag in der Theorie ganz toll klingen, praktisch wird da aber kein Nutzwert drin liegen. Du solltest Deine Filterordnung sinnvoll den Anforderungen entsprechend begrenzen. Der Ansatz "Viel hilft viel" ist hier nicht sinnvoll. Auch sehr hohe Ordnungen verhindern nicht, dass Du Dir bei jeder Faltung neben Deinem SI auch noch ein Rechteck mit in Dein Signal faltest. Selbst mit Ordnungen im Bereich von 10000 und mit double-Präzision wirst Du kaum mehr als 30 dB Sperrdämpfung erreichen. Vom Rechenaufwand pro Sample rede ich mal gar nciht erst. Fensterfunktionen oder andere Designmethoden sind hier eher das Mittel der Wahl.

PhilippHToner

Okay, das stimmt natürlich auch wieder.
Ich habe ja jetzt das, was ich wollte. Ich fenstere sinc. Je nach Belieben kann man die Breite des Fensters anpassen. Diese wird beim der Faltung berücksichtigt, sodass keine unnötigen Operationen ausgeführt werden. Höhere Genauigkeiten sind aufgrund der vorhandenen Daten eh nicht möglich.

Die Optimierung daher:
Die Berechnungszeit von 1sec*192000Samples/s auf 1sec*44100Samples/s ohne Optimierung(keinere Fensterung, keine Threads, kein Sinc-Table) würde auf meinem Rechner pro 1 Sample von den 44,1k Samples ~2sec brauchen. Das wären 44100*2 sec = 1 Tag.

Deshalb implementiere ist das massivparallel mit OpenCL. Neben der diskreten Implementierung werde ich das Ganze auch für die spektrale Seite implementieren, denn eine Faltung im Ortsbereich ist eine Multiplikation im Frequenzbereich.

PhilippHToner

Tachyon schrieb:

Erstmal eine Frage: Ist das Umtastverhätnis konstant?

Also ich taste nicht wie in dem ersten Beispiel von x kHz nach x kHz ab, weil ich mir das sonst sparen könnte. Es sind die gängigen Raten aus der Audiowelt denkbar:
48kHz Studioqualität => 44,1kHz Audio-CD.

Tachyon

PhilippHToner schrieb:

Okay, das stimmt natürlich auch wieder.
Ich habe ja jetzt das, was ich wollte. Ich fenstere sinc. Je nach Belieben kann man die Breite des Fensters anpassen. Diese wird beim der Faltung berücksichtigt, sodass keine unnötigen Operationen ausgeführt werden. Höhere Genauigkeiten sind aufgrund der vorhandenen Daten eh nicht möglich.

Die Optimierung daher:
Die Berechnungszeit von 1sec*192000Samples/s auf 1sec*44100Samples/s ohne Optimierung(keinere Fensterung, keine Threads, kein Sinc-Table) würde auf meinem Rechner pro 1 Sample von den 44,1k Samples ~2sec brauchen. Das wären 44100*2 sec = 1 Tag.

Deshalb implementiere ist das massivparallel mit OpenCL. Neben der diskreten Implementierung werde ich das Ganze auch für die spektrale Seite implementieren, denn eine Faltung im Ortsbereich ist eine Multiplikation im Frequenzbereich.

Hmm, eigentlich ist sowas ohne Optimierungen in vielfacher Echtzeit möglich.
Ich arbeite hier mit einem adaptiven Umtastfilter (das Umtastverhätnis ist Zeitvariant), und das Ding braucht um von 6kHz-50kHz variabel auf 44,1 kHz Fix zu kommen auf relativ alter Hardware (Pentium M, 1.8 GHz) ca. 3% CPU. Und Artefakte sind hier selbst mit schwachsinnig hoch aufgelösten Spektralschätzern nicht zu finden.

Was Du da tust ist aus meiner Sicht ist völlig unpraktikabel. Was willst Du erreichen?

PhilippHToner

Ich habe schon etliche Mal von Tonstudiogurus gehört, dass man die Finger von den billigen Echtzeit-Resamplern lassen soll und ich kann akustisch bestätigen, dass der Resampler von Cubase 5 von 48kHz auf 44,1kHz Artefakte produziert und diese hörbar sind. Zu jedem Bassdrumkick in dem Song bekomme ich eine hohe Pfeiffrequenz, die ich im Studio ohne resampling nicht bekomme.

Jetzt ist die Frage: Ist der Resampler billig-Software oder ist echtzeit-resampling einfach viel ungenauer als post-resampling?

hustbaer

Mal ne doofe Frage nebenbei: macht es überhaupt Sinn Audio-Daten mit einem sinc Kernel zu resampeln?
Das müsste doch zu massivem Ringing führen.

Ich hatte eigentlich den Eindruck, dass bei Audio eher Filter mit weniger steiler Flanke verwendet werden.

krümelkacker

PhilippHToner schrieb:

Ich hatte halt gedacht, wenn der sin(a) zuerst zu fehleranfällig wird, rechne ich halt mit sin( fmod(a, 2*M_PI)). Somit ist die math. Reihe doch genauer?

Wenn, dann bitte so: sin(fmod(x,2)*M_PI)

Mit steigendem x wird aber trotzdem die Genauigkeit von fmod(x,2) zunehmend schlechter; denn der absolute Fehler von x bleibt bei fmod(x,2) ungefähr erhalten.

Ich verstehe, dass es Dir um eine "möglichst fehlerfreie" Bandbeschränkung geht. Das Wort "fehlerfrei" ist aber schon unglücklich gewählt in diesem Kontext. Die Länge des "optimalen" Faltungskerns ist ein Zielkonflikt zwischen

steiler Flanke im Spektrum
kurzen "Ringing" in der Zeitdomäne.
hoher "stopband rejection" (Stärke der Dämpfung im Frequenzband, was man rausfiltern will)

In der Praxis sollte man einfach genug "headroom" im begrenzten Spektrum haben. Das heißt, man wählt eine Abtastrate, die meinetwegen 2,5 mal so hoch ist wie die höchste Frequenz die einen noch interessiert (statt 2,0+epsilon). Und dieser "headroom" erlaubt es dir, Tiefpassfilter mit kurzer Impulsantwort und dementsprechend einer flacheren Flanke zu nutzen.

Wenn man z.B. ein digitales Audiosignal von 96 kHz zu 44.1 kHz konvertiert, würde man einen Tiefpassfilter mit einem Passband bis 20 kHz und einer Flanke zwische 20-22.05 kHz nehmen. Länge und Art des optimalen Fensters hängt ein bißchen davon ab, was Du genau machen willst. Im Audiobereich sollte man bei Filtern mit einer Sperrdämpfung von mit mindestens 100 dB arbeiten. => Kaiser-Fenster mit entsprechendem Parameter.

Rezept zur Fensterfindung:

sich überlegen, was für eine "stopband rejection" man braucht und sich dementsprechend ein Fenster raussuchen, was dies bietet.
entsprechende Fensterlänge in Abhängigkeit der gewünschten Flankensteilheit wählen, Fensterlänge und Flankenbreite sind umgekehrt proportional zueinander

Super finde ich die Übersicht bei Wikipedia:
http://en.wikipedia.org/wiki/Window_function

Edit: Wieder etwas gelernt: stopband rejection nennt man im Deutschen Sperrdämpfung. Danke Tachyon!

Tachyon

PhilippHToner schrieb:

Ich habe schon etliche Mal von Tonstudiogurus gehört, dass man die Finger von den billigen Echtzeit-Resamplern lassen soll und ich kann akustisch bestätigen, dass der Resampler von Cubase 5 von 48kHz auf 44,1kHz Artefakte produziert und diese hörbar sind. Zu jedem Bassdrumkick in dem Song bekomme ich eine hohe Pfeiffrequenz, die ich im Studio ohne resampling nicht bekomme.

Jetzt ist die Frage: Ist der Resampler billig-Software oder ist echtzeit-resampling einfach viel ungenauer als post-resampling?

Das ist keine Billigsoftware, sondern ein selbstgeschriebener Algorithmus für HiFi-Audiostreaming für eine Kundengruppe die hier recht hohe Anforderungen hat. Preislich dürfte die Zielsysteme deutlich über professionellen Tontechniksystemen liegen.

Tachyon

krümelkacker schrieb:

Wenn man z.B. ein digitales Audiosignal von 96 kHz zu 44.1 kHz konvertiert, würde man einen Tiefpassfilter mit einem Passband bis 20 kHz und einer Flanke zwische 20-22.05 kHz nehmen. Länge und Art des optimalen Fensters hängt ein bißchen davon ab, was Du genau machen willst. Im Audiobereich sollte man bei Filtern mit einem Stopband mit mindestens 100 dB arbeiten. => Kaiser-Fenster mit entsprechendem Parameter.

Jo, das ist sinnvoll. Wobei sich die Fensterlänge aus dem gewünschten Transitionsband und der gewünschten Sperrdämpfung ergibt.

krümelkacker

heheh, genau. Habe mein Vorgänger Post noch um ein "Rezept" ergänzt. Ich bin nur nicht mit den deutschen Fachwörtern vertraut. Ich nehme mal an, dass Sperrdämpfung das ist, was ich vorhin "stopband rejection" nannte.

Tachyon

krümelkacker schrieb:

heheh, genau. Habe mein Vorgänger Post noch um ein "Rezept" ergänzt. Ich bin nur nicht mit den deutschen Fachwörtern vertraut. Ich nehme mal an, dass Sperrdämpfung im Englischen "stopband rejection" heißt.

Das kann man 1:1 mit "stopband attenuation" übersetzen.

Tachyon

Für die von krümelkacker genannten Randbedingungen brauchst Du mit der Fenstermethode und einem Kaiser-Fenster übrigens eine Ordnung von ca. 140. Damit solltest Du dann auch "etwas" schneller umtasten können als 1 Sample/2 s.