DE102005015456A1

DE102005015456A1 - Verfahren zur Bestimmung der zeitlichen Lage eines Wellenpakets sowie Flussmessgerät

Info

Publication number: DE102005015456A1
Application number: DE102005015456A
Authority: DE
Inventors: Eckard M. A. Glaser; Daniel Scheder
Original assignee: Viasys Healthcare GmbH
Current assignee: Vyaire Medical GmbH
Priority date: 2005-04-04
Filing date: 2005-04-04
Publication date: 2006-10-05
Also published as: EP1875380A2; EP1875380B1; WO2006105761A3; DE112006001498A5; US20090276167A1; US8103461B2; WO2006105761A2

Abstract

Die Erfindung betrifft ein Verfahren zur Bestimmung der zeitlichen Lage eines Wellenpakets (31) in einem Flussmessgerät sowie ein Flussmessgerät. Das Verfahren umfasst das Abtasten des Wellenpakets zu einer Vielzahl von Zeitpunkten, wobei ein Messwert zu jedem Zeitpunkt entsteht. Anschließend wird eine Summe von Produkten (43, 54, 55, 63, 64) berechnet, wobei jedes Produkt für einen bestimmten Zeitpunkt aus der Vielzahl von Zeitpunkten berechnet wird und jedes Produkt das Produkt aus einem Wert einer Vergleichsfunktion (44, 53, 56, 61, 62) zu dem bestimmten Zeitpunkt und dem Messwert zu dem bestimmten Zeitpunkt ist. Schließlich wird die zeitliche Lage des Wellenpakets aus der Summe von Produkten berechnet.

Description

Die Erfindung bezieht sich auf die Signalaufbereitung, um einen Laufzeitunterschied exakt zu messen. Die Erfindung bezieht sich ferner auf Flussmessgeräte, insbesondere Flussmessgeräte zum Erfassen des Luftflusses beim Ein- und Ausatmen, die den Mitnahmeeffekt von Ultraschall ausnutzen. Solche Flussmessgeräte werden auch als Spirometer bezeichnet.
Es sind eine Reihe von Geräten bekannt, die den Mitnahmeeffekt von Schallwellen in strömenden Medien, insbesondere Luft ausnutzen. Hierbei gibt es im wesentlichen zwei Bauformen, die in den 11 und 12 dargestellt sind.
Das in 11 dargestellte Flussmessgerät 1 umfasst ein Rohr 2 sowie einen ersten Ultraschallwandler 4 und einen zweiten Ultraschallwandler 5. Das Rohr hat einen Durchmesser d und eine Achse 3. Die beiden Ultraschallwandler 4 und 5 sind auf einer zweiten Achse 6 angeordnet, sodass jeder der beiden Wandler den vom jeweils anderen Wandler erzeugten Ultraschall auffängt. Die zweite Achse 6 schließt mit der Achse 3 des Rohres einen Winkel α ein. Bei ruhendem Medium im Rohr 2 gilt für die Schalllaufzeiten t₄₅ und t₅₄ von Ultraschallwandler 4 zu Ultraschallwandler 5 bzw. von Ultraschallwandler 5 zu Ultraschallwandler 4
wobei L der Abstand der beiden Ultraschallwandler und c die Schallgeschwindigkeit im ruhenden Medium bezeichnen.
Durch eine Strömung im Rohr 2 nach rechts verkürzt sich die Schalllaufzeit t₄₅ von Ultraschallwandler 4 zu Ultraschallwandler 5. Gleichzeitig verlängert sich die Schalllaufzeit t₅₄ von Ultraschallwandler 5 zu Ultraschallwandler 4. Es gilt für ein rechteckförmiges Strömungsprofil:
Hierbei bezeichnet v die Geschwindigkeit des strömenden Mediums nach rechts. Durch Einführung des Laufzeitunterschieds Δt und der mittleren Laufzeit t₀ Δt ≡ t54 – t45 (4)
erhält man für die Strömungsgeschwindigkeit v:
Aus (2) und (3) kann die Schallgeschwindigkeit im ruhenden Medium berechnet werden:
Durch Einsetzen von Formel (7) in (6) erhält man:
Hierbei ist die Konstante C₁ der Umrechnungsfaktor zwischen Strömungsgeschwindigkeit v und Volumenstrom V .. In die Konstante C, geht im wesentlichen die Querschnittfläche von Rohr 2 ein, wobei r den Radius von Rohr 2 bezeichnet. C1 = C2πr2 (9)
Die Konstante C₂ ist für ein rechteckförmiges Strömungsprofil 1. In realiter bildet sich im Rohr ein Geschwindigkeitsprofil aus, sodass die Geschwindigkeit des strömenden Mediums von der Position im Rohr, vom Volumenfluss und der Viskosität des strömenden Mediums abhängt. Für typische Volumenflüsse, das Medium Luft und Rohrdurchmesser von zwei bis drei Zentimetern bildet sich an der Rohrwand ein laminarer oder turbulenter Strömungsbereich und in Rohrmitte ein turbulenter Strömungsbereich aus. Je größer der Volumenfluss, desto kleiner der laminare Bereich und desto größer der turbulente Bereich. Dieses Strömungsprofil ist für Nichtlinearitäten in der Sensorkennlinie verantwortlich.

Aus den Laufzeiten t₄₅ und t₅₄ oder dem Laufzeitunterschied Δt und der mittleren Laufzeit t₀ kann der Fluss im Rohr berechnet werden. Dieses Messprinzip ist unter anderem in EP 0 051 293 A1 , EP 0 243 515 A1 , EP 0 597 060 B1 , CH 669 463 A5 , WO 00/26618 A1 und EP 0 713 080 A1 beschrieben. Je kleiner der Winkel α, desto größer ist der Laufzeitunterschied bei gleichem Fluss.

Das in 12 dargestellte bekannte Flussmessgerät 11 umfasst ein U-förmiges Rohr, das aus einem waagrechten Schenkel 12 und zwei senkrechten Schenkeln 17 und 18 besteht. Darüber hinaus umfasst das Flussmessgerät 11 Ultraschallwandler 14 und 15, die auf einer Achse 13 angeordnet sind. Die Achse 13 ist gleichzeitig Achse von waagrechtem Schenkel 12. Vorteilhaft an diesem Aufbau ist, dass die Flussrichtung parallel zur Achse der beiden Ultraschallwandler verläuft und somit der Winkel α 0° und damit cos(α) = 1 wird. Nachteilig ist der komplizierte Aufbau sowie Wirbel, die in den Übergangsbereich inzwischen den senkrechten und dem waagrechten Schenkel entstehen. Eine solche Bauformen ist beispielsweise in WO 90/05283 A1 oder EP 1 279 368 A2 beschrieben.

Unabhängig von der Bauformen ist bekannt (z. B. WO 90/05283, EP 1 279 368 A2 ), die Ultraschallwandler mit akustischen Impedanzwandlern zu versehen, um die Schallübertragung zwischen dem strömenden Medium und dem Ultraschallwandler und umgekehrt zu verbessern.

Der Zustand und die Zusammensetzung des zu messenden Gases ändert sich während der Atmung. Die Umgebungsbedingungen werden als Ambient Temperature Pressure (ATP) bezeichnet, die in Body Temperature Pressure Saturated (BTPS) umgerechnet werden können. Ausgeatmetes Gas hat eine Temperatur von etwa 34°C bis 37°C und ist vollständig mit Wasserdampf gesättigt, entspricht also BTPS.

Darüber hinaus ist die Rundsing-Technik (Englisch: sing-around technique) beispielsweise aus der EP 0 713 080 A1 bekannt. Bei der Rundsing-Technik werden wiederholt Schallimpulse von einem Ultraschallwandler zum anderen Ultraschallwandler gesendet. Sobald ein Impuls empfangen wird, wird der nächste Impuls gesendet. Nach 1 bis 200 Impulsen wird die Richtung umgekehrt. Bei diesem Verfahren wird die Zeitauflösung der Schalllaufzeit auf Kosten der Messrate erhöht. Dies ist insbesondere bei Medien mit hoher Schallgeschwindigkeit, insbesondere Flüssigkeiten, oder falls die Laufzeit aus technischen Gründen nur mit starkem Rauschen behaftet bestimmt werden kann, von Bedeutung.

Zur Analyse insbesondere von sich zeitlich verändernden Funktionen oder Signalen sind beispielsweise die Korrelationsanalyse, die Fourieranalyse oder – transformation, die Kurzzeitfouriertransformation sowie die Wavelettransformation bekannt.

Der Korrelationskoeffizient φ zweier Funktionen f(t) und g(t), die von der Zeit t abhängen und im Intervall von T₁ bis T₂ definiert sind, ist definiert als:

Der Korrelationskoeffizient kann als Maß für die Ähnlichkeit zweier Funktionen interpretiert werden.

Ferner ist die Kreuzkorrelationsfunktion Φ (KKF) definiert als zeitlicher Mittelwert des inneren Produkts der beiden um die Zeit τ gegeneinander verschobenen Funktionen f(t) und g(t):

Periodische Funktionen mit der Periodendauer (T₂ – T₁) können durch Fourier-Transformation vollständig als gewichtete Summe von Sinus- und Kosinusfunktionen dargestellt werden. Die Fourierkoeffizienten können als Korrelationsfaktor einer periodischen Funktionen und den Basisvektoren des Fourierraums aufgefasst werden.

Die Fourier-Transformation ist prädestiniert, globale Aussagen über ein Signal f(t) zu treffen, da sich ihre Basisfunktionen über den gesamten Zeitbereich erstrecken und keine zeitliche Lokalisierbarkeit zulassen.

Zur Analyse von transienten oder zeitlich begrenzten Signalen wurde die Kurzzeitfouriertransformation (STFT) entwickelt, die eine Entwicklung einer eindimensionalen Zeitfunktionen in die zweidimensionale Zeit-Frequenz-Ebene entspricht. Dabei wird das Signal in als quasistationär angenommene Bereiche aufgeteilt und diese werden unabhängig voneinander einer gefensterten Fourier-Transformation unterzogen (Jürgen Niederholz, Anwendungen der Wavelet-Transformation in Übertragungssystemen, Universität Duisburg, 1999 (http://www.ub.uni-duisburg.de/diss/diss0016/inhalt.htm). Die Fensterfunktionen im Zeitbereich bewirkt ebenfalls eine Fensterung im Frequenzbereich, so dass die STFT über den ganzen Frequenzbereich eine konstante effektive Bandbreite Δω und über den gesamten Zeitbereich eine effektive zeitliche Ausdehnung Δt besitzt.

Die konstante effektive Bandbreite und die Beschränkung auf Sinus- und Kosinusfunktionen wurde bei der STFT immer noch als Einschränkung empfunden. Beide Beschränkungen werden bei der Wavelettransformation aufgehoben. Eine allgemeine Einführung in die Wavelettransformation findet sich in (Antje Ohlhoff, Anwendung der Wavelettransformation in der Signalverarbeitung, Universität Bremen, 1996). Die Funktion Ψ(t) stellt ein gültiges Wavelet dar, wenn eine reelle, endliche Konstante C_Ψ existiert, wobei Ψ ~(ω) die Fouriertransformierte der Funktionen Ψ(t) ist.

Eine wesentliche Idee der Wavelettransformation besteht darin, die Funktion Ψ(t) schon geschickt zu wählen und dem zu lösenden Problem anzupassen, so dass die zu analysierenden Funktionen mit einer geringen Anzahl von Wavelets beschrieben werden können. Anders ausgedrückt soll sich nur eine geringe Anzahl von Koeffizienten der Basisfunktionen nennenswert von 0 unterscheiden.

Ausgehend von einem Basiswavelet oder Motherwavelet werden die übrigen Basisfunktionen, die eine vollständige Basis für den Funktionsraum bilden, mit einem Translationsparameter b und einem Dilatationsparameter a gemäß Gleichung (13) erzeugt.

Durch die Gleichung (13) wird die so genannte Güte

über den Frequenzbereich konstant gehalten, der im Kontext der Wavelettransformation auch als Skalierungsbereich bezeichnet wird. Hierbei ist ω_M eine Mittenfrequenz. In der Natur treten viele Phänomene auf, bei denen die Güte konstant ist.

Es ist Aufgabe der Erfindung durch ein intelligentes Auswerteverfahren den Einsatz weniger teurer Bauelemente zu ermöglichen und/oder die Genauigkeit der Laufzeitmessung zu erhöhen.

Diese Aufgabe wird durch die Lehre der unabhängigen Ansprüche gelöst.

Bevorzugte Ausführungsformen der Erfindung sind Gegenstand der Unteransprüche.

Vorteilhaft an der Verwendung einer gewichteten Summe von Messwerten eines Wellenpakets zur Bestimmung der zeitlichen Lage des Wellenpakets ist, dass hierdurch ein großer Teil der Energie des Wellenpakets ausgewertet wird, womit die zeitliche Unschärfe klein wird und somit die zeitliche Lage sehr genau bestimmt werden kann.

Vorteilhafterweise wird zum Gewichten eine Vergleichsfunktion verwendet, die der Form des Wellenpakets ähnelt. Eine Sinus- und Kosinusfunktion ähnelt zumindest einer Periode oder wenigen Perioden des Ultraschallpakets.

Die Verwendung einer Sinus- und einer Kosinusfunktion als Vergleichsfunktionen ist vorteilhaft, weil hierdurch die Lage des Wellenpakets innerhalb einer Periode mit Hilfe einer Arcustangensfunktion bestimmt werden kann, ohne weitere Vergleichsfunktionen zu berechnen.

Die Verwendung einer exponentiellen Hüllkurve, die einer Gauß'schen Glockenkurve ähnelt, passt die Gesamtform der Vergleichsfunktion an die Form eines Wellenpakets an. So ist es möglich, praktisch die gesamte Energie des Wellenpakets auszuwerten und die zeitliche Lage mit einer optimalen Genauigkeit zu bestimmen.

Die iterative Berechnung der zeitlichen Lage eines Wellenpakets reduziert die Anzahl der notwendigen Berechnungen von Produktsummen und kann so Rechenzeit sparen.

Um zu prüfen, ob die zeitliche Lage des Wellenpakets nicht um eine oder mehrere Perioden verschoben ist, können die Vergleichsfunktionen um eine oder mehrere Perioden verschoben werden und an diesen Stellen Produktsummen berechnet werden. Zu diesem Zweck kann auch die geometrische Summe zweier Summenprodukten berechnet werden, was mit einem noch geringeren Rechenaufwand verbunden ist. Zum gleichen Zweck können auch die Messwerte direkt mit Bruchteilen von Maxima und/oder Minima verglichen werden.

Eine Maximumsuche in Produktsummen von Messwerten mit einer unterschiedlich weit verschobenen Vergleichsfunktion erlaubt es, den Programmcode kurz zu halten.

Eine quadratische Interpolation mittels einer Parabel ist gut geeignet zur Maximumsuche. Außerdem stimmt der Verlauf einer Parabel in der Nähe des Scheitels gut mit dem Verlauf eine Sinus- oder Kosinusfunktion bei einem Maximum überein, da die Taylorentwicklung zweiter Ordnung hier ausreichend ist.

Im Folgenden werden bevorzugte Ausführungsformen der Erfindung unter Bezugnahme auf die beiliegenden Zeichnungen näher erläutert. Dabei zeigen:
1 ein Flussdiagramm eines TiDO Algorithmus;
2 ein Flussdiagramm eines TriTraX Algorithmus;
3 ein Flussdiagramm eines Teiles eines WaveTraX Algorithmus;
4 ein Flussdiagramm eines zweiten Teiles des WaveTraX Algorithmus;
5 ein typisches empfangenes Ultraschallpaket beim TiDO, TriTraX oder WaveTraX Algorithmus;
6 die Anregungen beim P-SEX Algorithmus;
7 ein typisches empfangenes Ultraschallpaket beim P-SEX Algorithmus;
8 eine Referenzfunktion beim P-SEX Algorithmus;
9 eine Kreuzkorrelationsfunktion zwischen einem Ultraschallpaket und einer Referenzfunktion beim P-SEX Algorithmus;
10 eine Bauform eines bekannten Spirometers mit einer erfindungsgemäßen Beschattung; und
11 eine Bauform eines bekannten Spirometers; und
12 eine weitere Bauform eines bekannten Spirometers.
Wie sich aus obigen Formeln, insbesondere Formel (8) ergibt, muss insbesondere der Laufzeitunterschied Δt und die mittlere Laufzeit t₀ möglichst exakt gemessen werden. Die anderen Größen, die in den Volumenstrom V eingehen, also C₁, C₂, L, α, und/oder r, sind weitgehend konstant und können entweder direkt oder über Kalibrationsmessungen bestimmt werden. Auch die Nichtlinearitäten, die dazu führen, dass C₁ und C₂ von V und damit von Δt abhängen, kann man mit Eichmessungen und Anpassung von Modellkurven wie beispielsweise Polynomen an die Messwerte in den Griff bekommen.
Es verbleibt also das technische Problem, den Laufzeitunterschied Δt oder die Laufzeiten t₄₅ und t₅₄ ausreichend genau zu messen. Dieses Problem ist äquivalent dazu, die zeitliche Lage von Ultraschallpaketen zu bestimmen, die ja Anfang und Ende von Δt festlegen.
Bei den hier beschriebenen Ausführungsformen werden Ultraschallwandler SC.131.XLV (SECO Sensor Consult GmbH) eingesetzt. Für diese Ultraschallwandler wird eine typische Mittenfrequenz von 310kHz, eine typische Bandbreite (±3dB) von ±40kHz sowie eine Schallkeulenbreite von 11° (–3dB) angegeben. Es wird Ultraschall der Frequenz f_US=312,5kHz verwendet, weil es vorteilhaft aber nicht notwendig ist, wenn die Abtastrate f_S ein ganzzahliges Vielfaches der Ultraschallfrequenz f_US ist Ultraschall unterliegt bei Normalbedingungen praktisch keiner Dispersion. Für die angestrebte Genauigkeit der Flussmessung muss die zeitliche Lage eines Ultraschallpakets auf etwa eine Nanosekunde genau bestimmt werden. Es ist mit hohen Kosten verbunden, das von einem Ultraschallwandler gelieferte Signal mit einer Abtastrate von einem Gigahertz abtasten und dann das Maximum in den Abtastwerten suchen. Trotzdem wird diese Ausführungsform als eine erfindungsgemäße Ausführungsform betrachtet.
Es hat sich gezeigt, dass die im folgenden beschriebenen Algorithmen die geforderte zeitliche Genauigkeit von 1ns bereits bei Abtastraten f_S = 5 MHz liefern. Mit moderatem technischen und finanziellen Aufwand sind derzeit Abtastraten bis 25 MHz erreichbar.
Das vom Analog-Digital-Wandler 94 oder 95 gelieferte Signal enthält Messwerte W_j die einen zeitlichen Abstand von 1/f_S aufweisen. Ein typisches Ultraschallpaket ist in 5 gezeigt. Die lokalen Maxima sind mit M1 bis M8 durchnummeriert. M5 markiert das globale Maximum.
Die Messwerte W_j sind rauschbehaftet und weisen einen unbekannten Offset auf. Deshalb liefert die Suche nach lokalen Maxima oder Nullstellen ohne vorherige Signalaufbereitung keine brauchbaren Ergebnisse, selbst wenn die Abtastrate ausreichend hoch wäre.
TiDO Algorithmus
Der Time Domain Only (TiDO) Algorithmus beruht im wesentlichen auf einer Bestimmung der Phasenlage einer einzelnen Schwingungsperiode im Wellenpaket. Die meiste Energie und damit den größten Rauschabstand (SNR, signal to noise ratio) weist die Schwingung mit der größten Amplitude auf, wobei das additive Rauschen als konstant über das gesamte Wellenpaket angenommen wird. Die Laufzeit wird in mehreren Schritten bestimmt:
1. Im ersten Schritt 41 wird das absolute Maximum M = W_jM an der Stelle j_M bestimmt. Der erste Schritt wird nur dann durchgeführt, wenn beispielsweise nach dem Einschalten des Geräts kein vorausgehender Messwert τ_M,vorig zur Verfügung steht. Falls ein vorausgehender Messwert zur Verfügung steht, kann der Schritt 41 durch den Schritt 48 ersetzt werden, der weiter unten beschrieben wird.
Die maximal zu erwartete Änderung der Laufzeiten t₄₅ und t₅₄ liegt im Bereich von ±20μs, wobei die Laufzeiten etwa 160μs betragen. Einschließlich der Länge des Schallpakets von etwa 30μs müssen also in Schritt 41 etwa 70μs*f_S Datenpunkte W_i analysiert werden.
Unter Berücksichtigung der Periodizität des Schallpakets kann der Suchalgorithmus in Schritt 41 optimiert werden, so dass nicht alle Messpunkte überprüft werden müssen. Man beginnt bei einem Schätzwert für die Stelle des globalen Maximums j_M. Dies kann ähnlich zu Schritt 48 das in der vorausgehenden Messung bestimmte Maximum τ_m,vorig oder nach dem Einschalten des Geräts auch eine beliebige Stelle sein. Man initialisiert M mit W_m,vorig. Man prüft die Stelle j_M+1 rechts davon. Falls W_jM < W_jM+1 geht man in dieser Richtung weiter, bis W_jE < M * 0,99. Dabei wird das neue Maximum M gespeichert. Falls jetzt W_jM> 0,99M prüft man noch die andere Richtung nach links bis W_jA < 0,99M. Falls es ein neues Maximum gibt, wird das in M gespeichert. Als Zwischenergebnis erhält man so ein lokales Maximum M = W_jM. Im ungünstigsten Fall wurden bis jetzt etwas mehr als S/2 Stellen verglichen. Jetzt prüft man W_jM–nS und W_jM+nS, n = 1, 2, .... Das Maximum dieser Werte ist ein ausreichend guter Schätzwert für das absolute Maximum. Aufgrund der glockenförmigen Hüllkurve stellt der Schätzwert nicht notwendigerweise das Maximum dar, und es ist zu erwarten, dass der Schätzwert um so schlechter ist, je höher n ist.
2. Die Datenpunkte W_j im Intervall j∊[jM – S/2;jM + S/2[ (14)entsprechen nun der Schwingungsperiode mit der größten Amplitude, die mittels einer rechteckförmigen Fensterfunktionen in Schritt 42 zur weiteren Signalverarbeitung "ausgeschnitten" wird. Dabei ist S die Anzahl der Messwerte in einer Schwingungsperiode des Ultraschallpakets, wobei S = fS/fUS (15)
Für eine Abtastfrequenz f_S von 5MHz und f_US von 312,5kHz hat S einen Wert von 16.
3. Nun wird die zeitdiskrete zyklische Variante Y_τ der Kreuzkorrelationsfunktion in Schritt 43 berechnet:
Die Berechnung von Gleichung (17) erfolgt in Schritt 44. Die Werte f_i können offline, zum Beispiel während der Kompilierung oder Assemblierung berechnet werden und in einem Feld, das auch als Lookup-Tabelle bezeichnet wird, abgelegt werden. Aus diesem Feld werden die Werte f_i zur Laufzeit ausgelesen. Bei Verwendung von Formeln (16) und (17) müssen in der Lookup-Tabelle etwa zwei Perioden das Kosinus gespeichert werden. Durch explizite Rückführung auf die Grundperiode gemäß Gleichung (17a) kann die Lookup-Tabelle auf eine Periode beschränkt werden.
4. In Y_τ wird nur das Maximum an der Stelle τ_M in Schritt 45 gesucht. Es ist zu erwarten, dass τ_M ≈ j_M ist. Eine leichte Verschiebung kann sich durch Rauschen in den Messpunkten W_j ergeben. Die Schrittweite von Y_τ entspricht der von W_j und damit 1/f_S. Sie genügt damit noch nicht der erforderlichen Genauigkeit von 1ns.
5. Eine Möglichkeit zur genaueren Maximumsbestimmung ist, eine Parabel P_τ gemäß Formel (18) in Schritt 46 durch das Maximum Y_τM und die beiden Nachbarpunkte Y_τM–1 und Y_τM–1 ZU legen und die Position des Scheitelpunkts τ_Max gemäß Formel (19) in Schritt 47 zu bestimmen. Pτ = aYτ2 + bYτ + c (18)
Die zeitliche Lage des Ultraschallpakets ergibt sich nun aus Gleichung (20):
Der TiDO Algorithmus würde bereits eine ausreichende Genauigkeit liefern. Es ist jedoch nötig, viele Datenpunkte für den ersten Schritt 41 zu analysieren, da sichergestellt werden muss, dass, auch bei maximaler Änderungsrate der Schalllaufzeit vor allem auf Grund einer Änderung der Strömungsgeschwindigkeit v des Mediums, sich das in Schritt 42 ausgewählte Intervall in den zu analysierenden Messpunkten befindet. Deshalb wird noch Schritt 48 eingeführt:
6. Das in der vorausgehenden Messung bestimmte Maximum τ_M,vorig wird in der nachfolgenden Messung als das im ersten Schritt 41 gefundene Maximum j_M angenommen (Schritt 48), wodurch die Maximumsuche in Schritt 41 überflüssig wird. Es ist dabei nicht wichtig, dass das Maximum M = WEM in der Mitte des in Schritt 42 ausgewählten Intervalls liegt. Auch die beiden Nachbarperioden enthalten nicht viel weniger Energie als die Periode mit dem Maximum. Wichtig ist jedoch, dass das Maximum in dem in Schritt 42 ausgewählten Intervall liegt. Dies ist der Fall, wenn sich die Phase zwischen zwei Wiederholungen der Schritte 42 bis 49 nicht um mehr als eine halbe Periode ändert. Da beim Prototyp Messraten von 1 kHz erreicht wurden und aus den Laufzeiten von etwa 160μs eine maximale Messrate von 6,25kHz abgeleitet werden kann, stellt dies keine wirkliche Einschränkung dar.
Diese Überlegungen gelten jedoch nur für kleine Flüsse, bei denen Wirbelzonen ausreichend klein sind. An der oberen Messgrenze von 20l/s werden bei einem Rohrdurchmesser von etwa 2 cm Geschwindigkeiten von 65 m/s oder 240 km/h erreicht. Bei hohen Flüssen wurde beobachtet, dass die Laufzeit zwischen zwei aufeinander folgenden Messungen so stark variiert, dass sich die Referenzstelle – also j_M beim TiDO Algorithmus oder j₉₃ beim TriTraX Algorithmus – zwischen zwei aufeinander folgenden Messungen um mehr als eine Ultraschallperiode ≈ 3μs unterscheidet. Dies scheint auf Wirbel zurückzuführen zu sein. Diese Wirbel führen zusammen mit der Signalverwehung zu einer um 20dB kleineren Amplitude des Empfangsignals bei einem Fluss von 20 l/s.
Aufgrund der Variation der Laufzeit bei hohen Flüssen kann es sinnvoll sein, in Schritt 45 das Maximum Y_τM mit den benachbarten Maxima Y_τM±nS n = 1, 2, ... oder W_τM mit W_τM±nS zu vergleichen. Sollte sich hierbei herausstellen, dass Y_τM±nS bzw. W_τM±nS das absolute Maximum ist, sollten die Schritte 42 bis 45 mit τ_M±nS als Schätzwert für j_M wiederholt werden.
Aufgrund der Variation der Laufzeit bei hohen Flüssen kann es bei allen in dieser Schrift beschriebenen Ausführungsformen ferner sinnvoll sein, die Laufzeiten aus mehreren aufeinander folgenden Messungen zu mitteln. Die Anzahl der Bemittelten Messungen kann mit dem Fluss zunehmen.
Nachteilig ist, dass sich aufgrund des Temperaturverhaltens und der Flussabhängigkeit der Richtcharakteristik der Ultraschallwandler die in 5 dargestellte Wellenform des Ultraschallpakets ändern kann. Deshalb, aufgrund von Rauschen, und des geringen Höhenunterschieds von 1 % zwischen den Maxima M5 und M6 kann das Maximum M6 zum absoluten Maximum werden. Bei starker Verwirbelung kann sogar M4 oder M7 zum globalen Maximum anwachsen. Deshalb erweist sich der TiDO Algorithmus als wenig robust und praxistauglich.
Nachteilig ist ferner, dass die Maxima der angefitteten Kosinusfunktion fi und der Parabel P_τ umso weiter auseinander liegen, je geringer die Abtastrate f_S ist.
TriTraX Algorithmus
Diese beiden Probleme soll der Trigonometric Transform by Cross (X) Correlation (TriTraX) Algorithmus lösen. Deshalb wird ein Punkt im Wellenpaket benutzt, der die Schwingungsperiode innerhalb des Ultraschallpakets zuverlässiger identifiziert. In 5 beträgt der Höhenunterschied zwischen den Maxima M4 und M5 etwa 14% der Höhe von M5 und ist zuverlässig detektierbar.
1. Im ersten Schritt 51 wird wieder das absolute Maximum M = Wj_M des Ultraschallpakets gesucht.
Aufgrund der Variation der Laufzeit bei hohen Flüssen kann es hierbei sinnvoll sein, in Schritt 51 das Maximum W_jM mit den benachbarten Maxima WτM±nS, n = 1, 2, ... zu vergleichen. Sollte sich hierbei herausstellen, dass W_τM±nS das absolute Maximum ist, wird j_M= j_M±nS gesetzt.
2. Nun wird der Punkt j₉₃ in bestimmt, dessen Ordinatenwert W_j93 am nächsten bei 0,93*M liegt oder diesen Wert zum ersten Mal in einem Wellenpaket überschreitet. Zur Optimierung des Signalrauschabstandes (SNR) wird die Periode mit der größten Amplitude für die weiteren Schritte im Algorithmus in Schritt 52 ausgewählt:
j₉₃ liegt etwa 1/8 Periode vor M5. Für die exakte Bestimmung der exakten, also interpolierten Position Maximums M5 ist es nicht wichtig, dass die ausgewählte Periode exakt von Minimum zu Minimum verläuft.
In anderen Ausführungsformen kann jedoch auch ein anderer Ordinatenwert als 93% des Maximums gewählt werden. Insbesondere kann dieser Ordinatenwert auch auf das Minimum, das arithmetische Mittel von Maximum und Minimum, die Differenz zwischen Maximum und Minimum und/oder den Gleichanteil des Wellenpakets bezogen sein. Ist der Ordinatenwert kleiner als der Gleichanteil des Wellenpakets, ist es besser, ein erstmaliges Unterschreiten des Ordinatenwerts als Bezugspunkt zu wählen.
3. Nun werden die Korrelationskoeffizienten φ_C und φ_S zwischen c_j und W_j bzw. s_j und W_j in Schritten 54 und 55 berechnet, wobei
Die Formeln (22) und (23) werden in Schritten 53 und 56 berechnet. Dies kann offline erfolgen, sodass die Werte während des Betriebs aus Feldern (Lookup-Tabellen) gelesen werden.
Die Berechnung der Korrelationskoeffizienten entspricht mathematisch der Bestimmung eines komplexen Koeffizienten φ_C + iφ_S einer diskreten Fourierentwicklung:
4. An Stelle dieser komplexen Darstellung kann auch der Betrag um die Phase φ der Schwingung angegeben werden, wobei die Phase φ mittels der arctan-Funktion berechnet werden kann. Diese ist jedoch nur in einem Intervall der Länge π eindeutig. Durch Berücksichtigung der Vorzeichen von φ_C und φ_S kann dieses Intervall auf 2π ausgedehnt werden: φ = arctan 2(φS;φC) (25)
In entsprechender Weise ist arc cot2(φ_C;φ_S) definiert.
Die Phase φ wird in Schritt 57 berechnet.
5. Die Laufzeit wird aus Formel 26 berechnet:
6. Analog zum oben beschriebenen TiDO Algorithmus wird wieder die Annahme in
Schritt 59 getroffen, dass
die Stelle des absoluten Maximums M = W_jM ist, um Schritt 51 aus Gründen der Rechengeschwindigkeit zu umgehen.
Wie beim TiDO Algorithmus kann zunächst ein lokales Maximum gesucht werden und dann im Abstand von ±nS benachbarte lokale Maxima überprüft werden, um das globale Maximum zu finden, wobei n eine ganze Zahl ist.
Sowohl der TiDO als auch der TriTraX Algorithmus haben den Nachteil, dass nur eine Periode des gesamten Ultraschallpakets ausgewertet wird, was zu einer suboptimalen Genauigkeit der Bestimmung der Laufzeit führt.
WaveTraX Algorithmus
Der Wavelet Transform by Cross (X) Correlation Algorithmus ermöglicht es, eine dem Signal angepasste Basis zu verwenden. Das Motherwavelet wird so gewählt, dass es möglichst genau einem Ultraschallpaket entspricht.
Aus einem nichtlinearen Fit über ein typisches Ultraschallpaket ergab sich das Mothersinus- und -cosinuswavelet zu:
Die Motherwavelets werden in den Schritten 61 und 62 berechnet. Die Ergebnisse können in Feldern abgelegt sein.
1. Zunächst wird in Schritt 63 die Lage des Wellenpaket mit Hilfe einer Kreuzkorrelationsfunktionen aus den Messwerten W_j mit dem Mothersinuswavelet ermittelt.
Die Kreuzkorrelationsfunktion K_WS,τ Erreicht ihren Maximalwert, wenn die Lage des Mothersinuswavelets und des abgetasteten Ultraschallpakets möglichst guter Übereinstimmung ist. Eine Maximumsuche in K_WS,τ liefert die Stelle τ_M des Maximums der Kreuzkorrelationsfunktion. Das Maximum der Kreuzkorrelationsfunktion ist gleichzeitig der Korrelationsfaktor φ_WS zwischen den Messwerten und dem Mothersinuswavelet: φWS = KWS,τM (30)
In diesem Zusammenhang soll schon jetzt darauf hingewiesen werden, dass in einer verfeinerten Ausführungsform in die Maximumsuche eine Iteration eingebaut werden kann, gemäß der in Gleichung (30) bei einem beliebigen r, also nicht unbedingt bei τ_M aber beispielsweise bei τ_M,vorig, dem τ_M der vorangehenden Messung, begonnen wird. Aus Gleichung (34) wird dann in Schritt 66 ein verbesserter Schätzwert τ_B für τ_M ermittelt. Ergebnis dieser Iteration ist ein lokales Maximum. Die Maximumsuche in den Schritten 72 bis 80 beschränkt sich dann auf lokale Maxima, deren Abstand zum gerade bestimmten Maximum sich im Wesentlichen aus der Frequenz des Ultraschalls ergibt.
2. Die exakte Phasenlage des Wellenpaket wird ähnlich dem TriTraX Algorithmus aus den Korrelationsfaktoren zwischen den Messwerten und dem Mothersinus- und -cosinuswavelet in Schritt 65 ermittelt. Dies ist äquivalent zur Ermittlung der Phasenlage aus der Phase eines komplexen Korrelationsfaktors. Der Korrelationsfaktor zwischen den Messwerten und dem Motherkosinuswavelet wird aus
in Schritt 64 bestimmt.
3. Aus den Korrelationsfaktoren wird die Phase φ_W bestimmt: φW = arc cot 2(φWC; φWS ) (32)
Die Laufzeit ergibt sich aus:
Steht bereits ein Schätzwert τ_S für den Index τ_M des Maximums zur Verfügung, so kann für die gleichen Messwerte W_j ein verbesserter Schätzwert τ_B aus
in Schritt 66 berechnet werden. φ_W ergibt sich aus Gleichung (32), wobei φ_WS und φ_WC aus Gleichungen (35) und (36) in Schritten 63 und 64 berechnet werden.
τ_B gibt nicht notwendigerweise exakt die Phasenlage des Ultraschallpakets an, weil die Exponentialfunktion in den beiden Motherwavelets eine Hüllkurve für die Sinus- bzw. Kosinus-Anteile darstellt und diese verzerrt. Stimmt die Hüllkurve des empfangenen Ultraschallpakets (vgl. 5) noch nicht mit dem Exponentialanteil der beiden Motherwavelets überein, weil der Abstand vom globalen Maximum |τ_S– τ_M| noch zu groß ist, ist nicht zu erwarten, dass der verbesserte Schätzwert τ_B exakt die Phasenlage des Ultraschallpakets angibt. Insofern kann der Betrag von Korrelationsfaktor φ_WC oder φ_W in Schritt 81 als Maß für die Konvergenz herangezogen werden.
Gemäß Gleichung (29) wird das Maximum des Korrelationsfaktors φ_WS mit dem Mothersinuswavelet (vgl. Gleichung (30)) bestimmt. Aufgrund der begrenzten Abtastrate weicht dieses Maximum optimalerweise um bis zu ±1/(2f_S) von der Phasenlage des Ultraschallpakets ab. Die tatsächlichen Abweichungen sind noch etwas größer. Die Nachkommastellen vom Index τ_M, der ja zunächst ganzzahlig ist, werden sozusagen mit Hilfe des Motherkosinuswavelets berechnet. Solange das Motherkosinuswavelet noch eine Verschiebung liefert, deren Betrag größer als ein Schwellenwert ist, kann man davon ausgehen, dass noch keine ausreichende Konvergenz erreicht wurde. Eine sinnvolle Wahl dieses Schwellenwerts ist 1/f_S. Für S = 16 ergibt sich eine betragsmäßige Obergrenze für φ_W von π/8 oder 22,5° und damit φWC/φWS > 0,414 (37)
Diese Obergrenze kann für den optimalen Fall für φ_W auf 11,25° und für φ_WC/φ_WS auf 0,199 abgesenkt werden.
Falls diese Obergrenze überschritten wird und deshalb ein weiterer Iterationsschritt durchgeführt werden soll, wird τ_B = τ_M zwischen den beiden Iterationsschritten in Schritt 68 auf einen ganzzahligen Wert gerundet. Mit τ_M wurde der interpolierten Index des Maximums von φ_WC bezeichnet. Gegen diesen konvergiert τ_B. Andernfalls müssten die Messwerte W_j interpoliert und die Funktionswerte der beiden Motherwavelets S_j und C_j neu berechnet werden. Insofern ist τ_S in Gleichung (34) ganzzahlig, wohingegen τ_B typischerweise Nachkommastellen aufweist.
In einer Ausführungsform, in der ein Mikroprozessor verwendet wird, in dem die Berechnung der Arkustangens- oder Arkuskotangensfunktion nicht effizient implementiert ist und deshalb lange dauert, können diese Funktionen in Schritt 65 mit einer geringen Genauigkeit berechnet werden. Eine hohe Genauigkeit ist erst erforderlich, wenn die Bedingung in Schritt 81 erfüllt ist, wenn also der ganzzahligen Anteil von τ_B= τ_M bekannt ist. Beim Prototypen beträgt S = 16, sodass in einen Quadranten lediglich 4 Schwellenwerte gelegt werden müssen, um den ganzzahligen Anteil von Index τ_B zu bestimmen. Dies zeigt, dass die Berechnung der Arkustangens- oder Arkuskotangensfunktion mit der in Schritt 65 erforderlichen Genauigkeit mittels einer Lookup-Tabelle implementiert werden kann.
Wie oben erwähnt, führt diese Iteration lediglich zu einem lokalen Maximum von φ_WS, das nicht notwendigerweise gleich dem globalen Maximum von φ_WS entspricht. Auch ist der Abstand zweier lokaler Maxima nicht immer exakt gleich der Periodendauer des Ultraschalls, weil die Exponentialfunktion die Sinus- und Kosinusfunktion verzerrt. Insofern kann die Maximumsuche dahingehend modifiziert werden, dass abwechselnd ein Iterationsschritt gemäß Gleichung (34) und eine Überprüfung der beiden benachbarten lokalen Maxima durchgeführt wird. Die beiden Korrelationskoeffizienten φ_WC und φ_WS wurden an der Stelle τ_S berechnet. Bei round(τ_B) werden die Korrelationskoeffizienten erst im nächsten Iterationsschritt berechnet. Die Funktionen round() rundet eine Fließkommazahl auf die nächste ganze Zahl. Da τ_S anfangs lediglich eine Überlappung zwischen dem Ultraschallpaket und den Wavelets sicherstellen muss, kann φ_WS am Anfang der Iterationen auch negativ sein. Insofern bietet sich an, zur Maximumsuche den Betrag B(τ) zu verwenden, dessen Berechnung beispielhaft in Schritt 67 dargestellt ist:
Ein Betrag wird aber auch in Schritten 73, 75, 77 und 79 berechnet, auch wenn dies an anderen Stellen erfolgt.
Die Zahl S, die einen gebrochenen Anteil aufweisen kann, gibt die Zahl von Messwerten in einer Ultraschallperiode an. Deshalb sollen die Beträge B(round(τ_S ± S)) in Schritten 73 und 77 berechnet und mit dem Betrag B(τ_S) verglichen werden. Falls B(round(τ_S ± S)) > B(τ_S), kann als Schätzwert gemäß den Schritten 76 bzw. 80 in Verbindung mit Schritt 68 für die nächste Iteration τ_S+ = round(τ_B ± S) verwendet werden.
In einer weiteren Ausführungsform können, wenn der Betrag B(round(τ_S ± S)) größer als der Betrag B(τ_S) ist, in dieser Richtung weitere Stellen round(τ_S ± πS) hinsichtlich des Betrags B(round(τs ± πS)) der Korrelationsfaktoren in Schritten 75 bzw. 79 untersucht werden. Diese Stellen liegen um ein ganzzahliges Vielfaches der Periodendauer des Ultraschalls von τ_S entfernt. Ziel der Schritte 75 und 76 ist es, möglichst bereits bei der nächsten Iteration beim globalen Maximum zu landen.
Solche Feinheiten sind bei einer hohen Messrate der Flussgeschwindigkeit für die bei menschlicher Atmung typischerweise maximal zu erwartenden Beschleunigungen im Gasstrom weniger wichtig. Dann unterscheiden sich die Phasenlagen zweier Ultraschallpakete zweier aufeinander folgender Messungen wenig. Deshalb liegt τ_M bereits ausreichend nah bei τ_M,vorig, der Stelle des Maximums bei der vorhergehenden Messung, so dass, wenn τ_M,vorig als erster Schätzwert verwendet wird, lediglich ein Iterationsschritt nötig ist, um τ_M zu bestimmen. In einem nicht ganz so optimalen Fall, wenn die Bedingung in Schritt 81 zunächst nicht erfüllt ist, sind zwei Iterationen notwendig, wobei die erste Iteration den ganzzahligen Anteil von τ_M bestimmt. Nachdem durch die erste Iteration die exponentiellen Hüllkurven des Ultraschallpakets und der beiden Motherwavelets optimal zur Deckung gebracht wurden, wird während der zweiten Iteration der gebrochene Anteil der Fließkommazahl τ_M exakt bestimmt. Ein Wechsel von einem lokalen zum globalen Maximum ist dann Idealerweise nicht notwendig, was in Schritt 71 verifiziert werden kann.
Um festzustellen, ob sich ein Schätzwert τ_S oder ein verbesserter Schätzwert τ_B sich bereits in der richtigen Ultraschallperiode befindet, kann der Betrag B(τ_S) in Schritt 71 mit einem Grenzwert C_B verglichen werden. B(τS) > CB (39)
Dieser Betrag muss über dem Grenzwert C_B liegen, der von der Flussgeschwindigkeit und damit vom Laufzeitunterschied Δt aber auch von anderen Parametern wie dem Alter der Ultraschallwandler abhängen kann. Aufgrund des Intensitätsabfalls bei hohen Flüssen wird C_B mit zunehmendem Fluss abnehmen.
Eine weitere Möglichkeit, bei dem WaveTraX Algorithmus die richtige Periode zu bestimmen, besteht in einer Kombination mit dem TriTraX Algorithmus. So kann zunächst der Index j₉₃ bestimmt werden, dessen Ordinatenwert W_j93 am nächsten bei 0,93*M liegt. Hierbei ist M das Maximum der Messdaten, die gerade ausgewertet werden (vgl. Schritt 51). Als Schätzwert für den Index j_M des Maximums kann round(τ_B) in Schritt 59 verwendet werden, wobei τ_B aus der vorangehenden Messung stammt. Nach Bestimmung von j₉₃ kann der erster Schätzwert τ_S aus
berechnet werden. Der WaveTraX Algorithmus läuft dann in Schritten 63 und 64 weiter.
P-SEX Algorithmus
Es hat sich auch der WaveTraX Algorithmus insbesondere bei einer wesentlichen Reduzierung der Signalamplitude und Signalverzerrungen aufgrund von Turbulenzen bei hohen Flüssen als nur unzureichend robust herausgestellt. Daher wurde ein neuer Algorithmus mit der Bezeichnung P-SEX (Phase-Shifted Excitation and Cross Correlation (X)) entwickelt.
Die Anregung der Ultraschallwandler erfolgt über eine in 6 dargestellte rechteckförmige Wechselspannung, deren Phase sich nach der Hälfte der Anregungszeit um 180° verschiebt oder invertiert wird, was äquivalent ist. Insgesamt umfasst die Wechselspannung acht Perioden. Die Frequenz der Wechselspannung f_US beträgt wieder 312,5kHz und entspricht damit etwa der Resonanzfrequenz der Ultraschallwandler.
Das bei ruhendem Medium empfangene Signal, das beispielhaft in 7 dargestellt ist, weißt ebenfalls einen Phasensprung auf. Man kann sich das empfangene Signal aus zwei zueinander phasenverschobenen Hälften vorstellen. Wird das Signal bei höheren Flüssen schlechter, so ist der Übergangsbereich der Phase nicht mehr wohl definiert.
1. Daher wird als Referenzfunktion für die Berechnung einer Kreuzkorrelationsfunktion nicht wie beim WaveTraX Algorithmus eine Funktion verwendet, die dem empfangenen Signal ähnelt. Vielmehr wird die in 8 dargestellte Referenzfunktion verwendet, bei der der Übergangsbereich „ausgeschnitten" ist. Es wurde durch Simulationen gefunden, dass eine in 8 dargestellte Referenzfunktion mit einer Periode cos, 2 Perioden Null, einer Periode –cos eine in 9 dargestellte Kreuzkorrelationsfunktion mit dem empfangenen Signal liefert, bei der das globale Maximum um etwa 23% über dem nächstniedrigerem Maximum liegt. Dadurch springt der P-SEX Algorithmus nicht mehr so leicht auf ein falsches Maximum. Bei der bei den Algorithmen TiDO, TriTraX und WaveTraX verwendeten Anregung von 5 Perioden, wodurch sich eine in 5 dargestellte, typische Signalform des empfangen Signals ergibt, beträgt der Unterschied gerade 1%.
2. Noch besser ist es, das linke der beiden betragsmäßig größten Minima, das im Folgenden als linkes Minimum τ_IM bezeichnet wird, zu betrachten. Der Unterschied zum betragsmäßig nächstkleineren Minimum beträgt fast 50% im Gegensatz zu den 23% bei Betrachtung des Maximums.
Es wurde aber beobachtet, dass das Minimum links vom linken Minimum betragsmäßig sehr groß wird – insbesondere größer als der Betrag des (bei geringen Flüssen globalen) Minimums rechts des linken Minimums τ_IM. Gleichzeitig führt die bei hohen Flüssen veränderte Form der P-SEX Kreuzkorrelationsfunktion dazu, dass das linke Minimum τ_IM zum globalen Minimum betragsmäßig anwächst und wohl definiert ist (9).
Es wird also in der Kreuzkorrelationsfunktion KKF gemäß Gleichung (41) das globale Minimum y₁ an der Stelle τ gesucht. Dann werden die beiden benachbarten Minima y₀ und y₂ gemäß Gleichungen (42) beziehungsweise (43) an den Stellen τ – S beziehungsweise τ + S betrachtet. Bei geringen Flüssen ist y₀ das linke Minimum. Bei hohen Flüssen ist y₁ das zum globale Minimum angewachsene linke Minimum und y₂ das bei geringen Flüssen globale Minimum. Die Auswahl des linken Minimums τ_IM erfolgt in Gleichung (44). y1 = Min(KKF) = KKF(τ) (41) y0 = KKF(τ – S) (42) y2 = KKF(τ + S) (43)
Im Ergebnis sucht dieser Algorithmus prinzipiell das linke Minimum. Der Laufzeitunterschied, ab dem das linke Minimum zum globalen Minimum angewachsen ist, beträgt etwa 5μs. Im Vergleich hierzu beträgt der Laufzeitunterschied bei einem maximalen Fluss von 20 l/s bei der aktuellen Dimensionierung des Rohrquerschnitts etwa 24μs. Damit lässt sich die grobe Position des Schallpakets gut bestimmen.
3. Zur exakten Phasenbestimmung wird in der Kreuzkorrelationsfunktion die cos-Periode (–π ... π), die auf das gesuchte Minimum folgt, ausgewählt: Wi = KKF(τM + i), i = 0, 1, 2, ..., S (45)
4. Die W_i werden analog zu Schritten 54 und 55 des TriTraX Algorithmus nach einer Kosinus- bzw. Sinus-Funktion entwickelt, wobei die Koeffizienten φ_C bzw. φ_S berechnet werden.
5. Anschließend wird entsprechend dem Schritt 57 in 2 und Gleichung (25) mit einer arctan2-Funktion die Phase φ berechnet. Dies ergibt eine Phasenbestimmung mit geringerem Fehler, als eine Periode aus dem Signal zu betrachten oder das Maximum der Kreuzkorrelationsfunktion durch eine Parabel anzufitten.
6. Schließlich wird die Laufzeit t aus Gleichung (46) berechnet:
Neben den oben beschriebenen Algorithmen sind weitere Abwandlungen denkbar. Insbesondere kann das empfangene Ultraschallsignal mit einer Referenzfunktion korreliert oder gefaltet werden, die die Inverse des Übertragungswegs darstellt. Je nach Definition des Übertragungswegs, also beispielsweise Leistungsverstärker-Ultraschallwandler-Luft-Ultraschallwandler oder Triggerimpuls-Rechteckgenerator-Leistungsverstärker-Ultraschallwandler-Luft-Ultraschallwandler kann das Ergebnis der Faltung oder Korrelation beispielsweise ein Impuls oder eine Rechteck-Wechselspannung sein. Besonders vorteilhaft sind auch Korrelationsfunktionen, die das empfangene Ultraschallsignal näherungsweise in eine Dirac'sche Deltafunktion verwandeln, wobei Idealerweise nur der Wert eines Zeitintervalls nennenswert von Null verschieden ist.
Zur Bestimmung des richtigen Intervalls kann eine andere Korrelationsfunktion eingesetzt werden als zur Bestimmung der Phase.
10 zeigt eine Bauform eines bekannten Spirometers aus 11 mit einer erfindungsgemäßen Beschattung. Herzstück der Steuerung ist der Mikroprozessor 9. Als Mikroprozessor kann insbesondere ein digitaler Signal Prozessor eingesetzt werden, der die Berechnung von Summenprodukten besonders schnell durchführen kann. Solche Summenprodukte kommen beispielsweise in Gleichungen (16), (29), (31), (35) und (36) vor. Der Mikroprozessor 9 steuert insbesondere Treiberschaltung 7, die die Wechselspannung für die beiden Ultraschallwandler 4 und 5 erzeugt. Da beide Wandler gleichzeitig senden sollen, ist lediglich eine Treiberschaltung erforderlich.
Insbesondere während des Empfangsbetriebs müssen die Wandler von der Treiberschaltung getrennt werden, was beispielsweise durch die Schalter 8 erfolgen kann, die aktiv durch den Mikroprozessor 9 ein- und ausgeschaltet werden. Diese Trennung kann auch dadurch erfolgen, dass der Ausgang der Treiberschaltung 7 hochohmig geschaltet wird. Eine dritte Möglichkeit besteht im Einsatz nichtlineare Bauelemente wie Dioden, wie er beispielsweise in EP 0 234 515 A1 , 3 offenbart ist.
Im Empfangsbetrieb werden die von den Ultraschallwandlern gelieferten elektrischen Signale durch die Analog-Digital-Wandler 94 und 95 digitalisiert und dem Mikroprozessor 9 zur weiteren Verarbeitung zugeführt. Bei einem Prototypen wurden DSPs der Baureihe TMS320C24x (Texas Instruments) eingesetzt, die einen 12-Bit ADC und eine 32-Bit CPU enthalten.
Insbesondere das Rohr 2 wird nicht als wesentlicher Teil der Erfindung angesehen, weil es selbst oder zumindest ein Einsatzrohr wegen der besseren Desinfizierbarkeit auswechselbar ausgeführt ist und deshalb getrennt vom Spirometer vertrieben wird. Außerdem kann die mittlere Flussgeschwindigkeit in Richtung der durch die Wandler 4 und 5 festgelegten Achse 6 ein sinnvolles Messsignal für ein Anemometer sein, ohne dass eine Fläche spezifiziert ist und folglich kein Volumenstrom bestimmt werden kann.
Obwohl die Lage der ausgesandten Ultraschallpakete leicht zu bestimmen ist, weil ja der Mikroprozessor 9 die Treiberschaltung 7 triggert, soll betont werden, dass das technische Problem darin besteht, den Laufzeitunterschied Δt, weniger die Laufzeiten t₄₅ und t₅₄ oder die mittlere Laufzeit to, exakt zu bestimmen.
Aufgrund der in der Zwischenzeit zur Verfügung stehenden Ultraschallwandler können kurze Ultraschallimpulse erzeugt werden. So beträgt die Länge des in 5 dargestellten Ultraschallpakets 20 bis 30μs. Das Abklingverhalten nach einem Ultraschallimpuls ist so gut, dass der gleiche Ultraschallwandler nach Verstreichen der Laufzeit eines Ultraschallpakets von etwa 160μs nach dem Aussenden eines Ultraschallimpulses für den Empfang eines Ultraschallpakets verwendet werden kann. Deshalb werden bei einer Ausführungsform gleichzeitig von beiden Ultraschallwandlern 4 und 5 oder 14 und 15 Ultraschallimpulse ausgesendet und etwa 160μs später von beiden Ultraschallwandlern wieder empfangen. So liegen die Messpunkte W_j für je eine Laufzeitmessung in beide Richtungen bereits nach gut 160μs vor, so dass man aus diesen Messdaten bereits den Laufzeitunterschied Δt und die mittlere Laufzeit to berechnen kann.
Soweit ausreichend Speicher und Rechenleistung zur Verfügung stehen, können sofort nach der Aufzeichnung von zwei ersten Ultraschallpaketen die beiden nächsten Ultraschallpakete gesendet werden. Während der Laufzeit der nächsten Ultraschallpakete kann die Auswertung der Aufzeichnung der ersten Ultraschallpakete erfolgen. Die Flussgeschwindigkeit, die aus den beiden ersten Ultraschallpaketen berechnet wird, kann dann kurz vor der Aufzeichnung der beiden nächsten Ultraschallpakete zur Verfügung stehen. So können die oben angegebenen Abtastraten von 6,25 kHz für die Geschwindigkeitsmessung erreicht werden.
Bei jeder der oben beschriebenen Ausführungsformen kann als Plausibilitätskontrolle der Unterschied der aus zwei aufeinander folgenden Messungen bestimmten Geschwindigkeit ausgewertet werden. Dieser muss betragsmäßig kleiner als eine maximale Änderungsrate der Geschwindigkeit multipliziert mit dem zeitlichen Abstand der beiden Messungen liegen. Falls dies nicht der Fall ist, kann eine Fehlermeldung ausgegeben werden und/oder der jüngere Messwert verworfen werden.
Obwohl die bevorzugte Ausführungsform oben im Zusammenhang mit Luft beschrieben wurde, kann die Erfindung bei jedem anderen Fluid eingesetzt werden, um eine Flussgeschwindigkeit oder ein Volumenstrom zu bestimmen.
Die Erfindung wurde zuvor anhand von bevorzugten Ausführungsformen näher erläutert. Für einen Fachmann ist jedoch offensichtlich, dass verschiedene Abwandlungen und Modifikationen gemacht werden können, ohne vom Geist der Erfindung abzuweichen. Deshalb wird der Schutzbereich durch die nachfolgenden Ansprüche und ihre Äquivalente festgelegt.
Mathematische Symbole:

α:

Winkel zwischen Achsen 3 und 6

c:

Schallgeschwindigkeit im ruhenden Medium

C_B, C₁, C₂:

Konstanten

S_j, C_j:

"Abtastwerte" des Mothersinus- und -cosinuswavelets

f_US:

Ultraschallfrequenz

f_S:

Abtastrate

f(t), g(t):

zeitabhängige Funktionen

j:

Index der Messwerte

j_M:

Index von M

L:

Abstand der Ultraschallwandler

M:

Maximum der W_j

φ_C, φ_S:

Korrelationskoeffizienten

φ_WS, φ_WC:

Korrelationskoeffizienten

φ_W:

Phase

Δt:

Laufzeitunterschied

t₀:

mittlere Laufzeit

t₄₅, t₅₄

Schalllaufzeiten

Φ:

Kreuzkorrelationsfunktion

r:

Radius von Rohr 2

S:

f_S/f_US (in einer Ausführungsform 16)

v:

Strömungsgeschwindigkeit des Mediums

V .:

Volumenstrom

W_j

Messwerte

Y_τ:

zeitdiskrete zyklische Variante

τ_M:

(interpolierter) Index

τ_S:

(ganzzahliger) Schätzwert für τ_M

τ_B:

verbesserter Schätzwert für τ_M

1, 11: Flussmessgerät
2: Rohr
3: Achse
4, 5, 14, 15: Ultraschallwandler
6: Achse
7: Treiberschaltung
8: Schalter
9: Mikroprozessor
94, 95: Analog-Digital-Wandler
12: waagrechten Schenkel
13: Achse
17,18: senkrechter Schenkel
31: Schallpaket
M1, M2, M3, M4, M5, M6, M7, M8: Maxima
41–83: Schritte

Claims

Verfahren zur Bestimmung der zeitlichen Lage eines Wellenpakets (31) in einem Flussmessgerät, insbesondere Spirometer mit: Abtasten (94, 95) des Wellenpakets zu einer Vielzahl von Zeitpunkten, wobei ein Messwert zu jedem Zeitpunkt entsteht; Berechnen einer Summe von Produkten (43, 54, 55, 63, 64), wobei jedes Produkt für einen bestimmten Zeitpunkt aus der Vielzahl von Zeitpunkten berechnet wird und jedes Produkt das Produkt aus einem Wert einer Vergleichsfunktion (44, 53, 56, 61, 62) zu dem bestimmten Zeitpunkt und dem Messwert zu dem bestimmten Zeitpunkt ist; und Berechnen (46, 47, 57, 58, 65–81) der zeitlichen Lage des Wellenpakets aus der Summe von Produkten.
Verfahren gemäß Anspruch 1, dadurch gekennzeichnet, dass die Vergleichsfunktion eine Sinusfunktion (56) ist, dass eine zweite Summe von zweiten Produkten berechnet wird (54), wobei jedes zweite Produkt für einen bestimmten Zeitpunkt aus der Vielzahl von Zeitpunkten berechnet wird und jedes zweite Produkt das Produkt aus einem Wert einer Kosinusfunktion (53) zu dem bestimmten Zeitpunkt und dem Messwert zu dem bestimmten Zeitpunkt ist, und dass die zeitliche Lage des Wellenpakets als Arcustangensfunktion (57) des Quotienten aus der Summe von Produkten dividiert durch die zweite Summe von zweiten Produkten berechnet wird.
Verfahren gemäß Anspruch 1, dadurch gekennzeichnet, dass die Vergleichsfunktion eine Sinusfunktion multipliziert mit einer Exponentialfunktion ist (61), wobei sich das Argument der Exponentialfunktion aus einer Potenz des mit einem Faktor gewichteten Betrags der Differenz des Zeitpunkts minus einem zeitlichen Offset ergibt.
Verfahren gemäß Anspruch 3, dadurch gekennzeichnet, dass eine zweite Summe von zweiten Produkten berechnet wird, wobei jedes zweite Produkt für einen bestimmten Zeitpunkt aus der Vielzahl von Zeitpunkten berechnet wird und jedes zweite Produkt das Produkt aus einem Wert einer zweiten Vergleichsfunktion zu einem bestimmten Zeitpunkt und dem Messwert zu dem bestimmten Zeitpunkt ist, wobei die zweite Vergleichsfunktion eine Kosinusfunktion multipliziert mit einer Exponentialfunktion ist (62), wobei sich das Argument der Exponentialfunktion aus der Potenz des mit dem Faktor gewichteten Betrags der Differenz des bestimmten Zeitpunkts minus dem zeitlichen Offset ergibt, und dass die zeitliche Lage des Wellenpakets als Arkuskotangensfunktion (65) des Quotienten aus der zweiten Summe von zweiten Produkten dividiert durch die Summe von Produkten berechnet wird.
Verfahren gemäß Anspruch 4, gekennzeichnet durch: Verschieben der Vergleichsfunktion sowie der zweiten Vergleichsfunktion (68) um den zeitlichen Abstand zwischen zwei Messwerten, so dass die Vergleichsfunktion möglichst gut mit der aus der Arkuskotangensfunktion berechneten zeitlichen Lage des Wellenpakets übereinstimmt; erneutes Berechnen (63, 64) der Summe von Produkten mit der verschobenen Vergleichsfunktion sowie der zweiten Summe von Produkten mit der verschobenen zweiten Vergleichsfunktion; und Berechnen eines besseren Schätzwerts für die zeitliche Lage des Wellenpakets als Arkuskotangensfunktion (65) des Quotienten aus der neu berechneten zweiten Summe von zweiten Produkten dividiert durch die Summe von Produkten. Verfahren gemäß Anspruch 4 oder 5, gekennzeichnet durch: Verschieben (73, 75, 77, 79) der Vergleichsfunktion um den zeitlichen Abstand zwischen zwei Messwerten, der etwa einem ganzzahligen Vielfachen einer Periode des Wellenpaktes entspricht, wobei die Vergleichsfunktion schon einmal verschoben wurde, um ihre Lage mit der zeitlichen Lage des Wellenpakets zur Deckung zu bringen; Berechnen der Summe von Produkten mit der verschobenen Vergleichsfunktion; Wiederholen des in diesem Anspruch definierten Verschiebens und Berechnens um eine Vielzahl von Summen von Produkten zu erhalten; und Auswählen der zeitlichen Verschiebung der Vergleichsfunktion, bei der die Summe von Produkten maximal wird.
Verfahren gemäß einem der Ansprüche 4 bis 6, dadurch gekennzeichnet, dass eine geometrische Summe (67) der Summe und der zweiten Summe berechnet wird und die geometrische Summe maximiert wird (73, 75, 77, 79), um die richtige Periode des Wellenpakets auszuwählen.
Verfahren gemäß einem der obigen Ansprüche, gekennzeichnet durch: Bestimmen des Messpunkts, an dem das Wellenpaket zum ersten Mal einen vorgegebenen Bruchteil des Maximums aller Messpunkte des Wellenpakets überschreitet. Auswählen (52) einer Periode des Wellenpakets um den bestimmten Messpunkt herum zur Berechnung der Summen von Produkten.
Verfahren gemäß einem der Ansprüche 1 bis 6, gekennzeichnet durch: Bestimmen des Messpunkts, an dem das Wellenpaket zum ersten Mal einen vorgegebenen Bruchteil des Minimums aller Messpunkte des Wellenpakets unterschreitet. Auswählen einer Periode des Wellenpakets um den bestimmten Messpunkt herum zur Berechnung der Summen von Produkten.
Verfahren gemäß einem der Ansprüche 1, 2, 3, 7, 8 oder 9, gekennzeichnet durch: Berechnen (43) einer Vielzahl von Summen von Produkten, wobei jeweils unterschiedliche Vergleichsfunktionen, aber die gleichen Messwerte verwendet werden, wobei die Vergleichsfunktionen um eine jeweils unterschiedliche Zeitdifferenz verschoben, sonst aber gleich sind; Bestimmen (45) des Maximums der Summen von Produkten; und Berechnen der zeitlichen Lage des Wellenpakets aus der Zeitdifferenz, um die die Vergleichsfunktion verschoben wurde, die zur Berechnung der maximalen Summe von Produkten verwendet wurde.
Verfahren gemäß Anspruch 10, gekennzeichnet durch: Legen (46) einer Parabel durch die maximale Summe von Produkten sowie die beiden benachbarten Summen, die aus Vergleichsfunktionen berechnet werden, die um den zeitlichen Abstand zweier benachbarter Messpunkte gegenüber der Vergleichsfunktion der maximalen Summe nach vorne und nach hinten verschoben sind; und Berechnen (47) der zeitlichen Lage des Wellenpakets aus der Lage des Scheitelpunkts der Parabel.
Verfahren gemäß einem der obigen Ansprüche, dadurch gekennzeichnet, dass der Ultraschallwandler, der das Wellenpaket erzeugt, mit einem Wechselspannungssignal angesteuert wird, das einen Phasensprung aufweist.
Verfahren gemäß Anspruch 12, dadurch gekennzeichnet, dass die Vergleichsfunktion aus einer Kosinus-Periode, zwei Perioden null, einer inversen Kosinus-Periode und sonst null besteht.
Flussmessgerät (1, 11) mit: einem ersten Schallwandler (4, 14); einem zweiten Schallwandler (5, 15), wobei die beiden Schallwandler aufeinander ausgerichtet sind, so dass der zweite Schallwandler (5, 15) einen Teil des vom ersten Schallwandler (4, 14) erzeugten Schalls aufnimmt und umgekehrt, wobei die beiden Schallwandler ein elektrisches Signal in Schall und umgekehrt wandeln; einer Treiberschaltung (7), die mit dem ersten Schallwandler (4, 5, 14, 15) verbunden ist, und den ersten Schallwandler zeitweise mit Wechselspannung zur Erzeugung von Schall versorgt; gekennzeichnet durch: einen Analog-Digital-Wandler (94, 95) zur Abtastung und Digitalisierung des vom zweiten Schallwandler ausgegebenen elektrischen Signals zu einer Vielzahl von Zeitpunkten, wobei ein Messwert zu jedem Zeitpunkt entsteht, wobei der Eingang des Analog-Digital-Wandlers (94, 95) mit dem zweiten Schallwandler (4, 5, 14, 15) verbunden ist, einen Mikroprozessor (9), der mit dem Ausgang des Analog-Digital-Wandlers (94, 95) verbunden ist und dem die Messwerte zugeführt werden, wobei der Mikroprozessor (9) die folgenden Schritte ausführt: Berechnen einer Summe von Produkten (43, 54, 55, 63, 64), wobei jedes Produkt für einen bestimmten Zeitpunkt aus der Vielzahl von Zeitpunkten berechnet wird und jedes Produkt das Produkt aus einem Wert einer Vergleichsfunktion (44, 53, 56, 61, 62) zu dem bestimmten Zeitpunkt und dem Messwert zu dem bestimmten Zeitpunkt ist; und Berechnen (46, 47, 57, 58, 65–81) der zeitlichen Lage des Wellenpakets aus der Summe von Produkten.
Flussmessgerät gemäß Anspruch 14, dadurch gekennzeichnet, dass die Vergleichsfunktion eine Sinusfunktion (56, 61) ist, dass der Mikroprozessor (9) eine zweite Summe von zweiten Produkten berechnet (54, 64), wobei jedes zweite Produkt für einen bestimmten Zeitpunkt aus der Vielzahl von Zeitpunkten berechnet wird und jedes zweite Produkt das Produkt aus einem Wert einer Kosinusfunktion (53, 62) zu dem bestimmten Zeitpunkt und dem Messwert zu dem bestimmten Zeitpunkt ist, und die zeitliche Lage des Wellenpakets als Arcustangensfunktion (57, 56) des Quotienten aus der Summe von Produkten dividiert durch die zweite Summe von zweiten Produkten berechnet.
Flussmessgerät gemäß Anspruch 14, dadurch gekennzeichnet, dass die Vergleichsfunktion eine Sinusfunktion multipliziert mit einer Exponentialfunktion ist (61), wobei sich das Argument der Exponentialfunktion aus einer Potenz des mit einem Faktor gewichteten Betrags der Differenz des Zeitpunkts minus einem zeitlichen Offset ergibt.
Flussmessgerät gemäß Anspruch 16, dadurch gekennzeichnet, dass der Mikroprozessor (9) eine zweite Summe von zweiten Produkten berechnet, wobei jedes zweite Produkt für einen bestimmten Zeitpunkt aus der Vielzahl von Zeitpunkten berechnet wird und jedes zweite Produkt das Produkt aus einem Wert einer zweiten Vergleichsfunktion zu einem bestimmten Zeitpunkt und dem Messwert zu dem bestimmten Zeitpunkt ist, wobei die zweite Vergleichsfunktion eine Kosinusfunktion multipliziert mit einer Exponentialfunktion ist (62), wobei sich das Argument der Exponentialfunktion aus der Potenz des mit dem Faktor gewichteten Betrags der Differenz des bestimmten Zeitpunkts minus dem zeitlichen Offset ergibt, und dass die zeitliche Lage des Wellenpakets als Arkustangensfunktion (65) des Quotienten aus der Summe von Produkten dividiert durch die zweite Summe von zweiten Produkten berechnet wird.
Flussmessgerät gemäß Anspruch 17, dadurch gekennzeichnet, dass der Mikroprozessor (9) ferner folgende Schritte ausführt: Verschieben der Vergleichsfunktion sowie der zweiten Vergleichsfunktion (68) um den zeitlichen Abstand zwischen zwei Messwerten, so dass die Vergleichsfunktion möglichst gut mit der aus der Arkuskotangensfunktion berechneten zeitlichen Lage des Wellenpakets übereinstimmt; erneutes Berechnen (63, 64) der Summe von Produkten mit der verschobenen Vergleichsfunktion sowie der zweiten Summe von Produkten mit der verschobenen zweiten Vergleichsfunktion; und Berechnen eines besseren Schätzwerts für die zeitliche Lage des Wellenpakets als Arkustangensfunktion (65) des Quotienten aus der neu berechneten Summe von Produkten dividiert durch die zweite Summe von zweiten Produkten.
Flussmessgerät gemäß Anspruch 17 oder 18, dadurch gekennzeichnet, dass der Mikroprozessor (9) ferner folgenden Schritte ausführt: Verschieben (73, 75, 77, 79) der Vergleichsfunktion um den zeitlichen Abstand zwischen zwei Messwerten, der etwa einem ganzzahligen Vielfachen einer Periode des Wellenpaktes entspricht, wobei die Vergleichsfunktion schon einmal verschoben wurde, um ihre Lage mit der zeitlichen Lage des Wellenpakets zur Deckung zu bringen; Berechnen der Summe von Produkten mit der verschobenen Vergleichsfunktion; Wiederholen des in diesem Anspruch definierten Verschiebens und Berechnens um eine Vielzahl von Summen von Produkten zu erhalten; und Auswählen der zeitlichen Verschiebung der Vergleichsfunktion, bei der die Summe von Produkten maximal wird.
Flussmessgerät gemäß einem der Ansprüche 17 bis 19, dadurch gekennzeichnet, dass der Mikroprozessor (9) eine geometrische Summe (67) der Summe und der zweiten Summe berechnet und die geometrische Summe maximiert (73, 75, 77, 79), um die richtige Periode des Wellenpakets auszuwählen.
Flussmessgerät gemäß einem der Ansprüche 14 bis 20, dadurch gekennzeichnet, dass der Mikroprozessor (9) ferner folgenden Schritte ausführt: Bestimmen des Messpunkts, an dem das Wellenpaket zum ersten Mal einen vorgegebenen Bruchteil des Maximums aller Messpunkte des Wellenpakets überschreitet. Auswählen (52) einer Periode des Wellenpakets um den bestimmten Messpunkt herum zur Berechnung der Summen von Produkten.
Flussmessgerät gemäß einem der Ansprüche 14 bis 20, dadurch gekennzeichnet, dass der Mikroprozessor (9) ferner folgenden Schritte ausführt: Bestimmen des Messpunkts, an dem das Wellenpaket zum ersten Mal einen vorgegebenen Bruchteil des Minimums aller Messpunkte des Wellenpakets unterschreitet. Auswählen einer Periode des Wellenpakets um den bestimmten Messpunkt herum zur Berechnung der Summen von Produkten.
Flussmessgerät gemäß einem der Ansprüche 14, 15, 16, 20, 21 oder 22, dadurch gekennzeichnet, dass der Mikroprozessor (9) ferner folgende Schritte ausführt: Berechnen (43) einer Vielzahl von Summen von Produkten, wobei jeweils unterschiedliche Vergleichsfunktionen, aber die gleichen Messwerte verwendet werden, wobei die Vergleichsfunktionen um eine jeweils unterschiedliche Zeitdifferenz verschoben, sonst aber gleich sind; Bestimmen (45) des Maximums der Summen von Produkten; und Berechnen der zeitlichen Lage des Wellenpakets aus der Zeitdifferenz, um die die Vergleichsfunktion verschoben wurde, die zur Berechnung der maximalen Summe von Produkten verwendet wurde.
Flussmessgerät gemäß Anspruch 23, dadurch gekennzeichnet, dass der Mikroprozessor (9) ferner folgenden Schritte ausführt: Legen (46) einer Parabel durch die maximale Summe von Produkten sowie die beiden benachbarten Summen, die aus Vergleichsfunktionen berechnet werden, die um den zeitlichen Abstand zweier benachbarter Messpunkte gegenüber der Vergleichsfunktion der maximalen Summe nach vorne und nach hinten verschoben sind; und Berechnen (47) der zeitlichen Lage des Wellenpakets aus der Lage des Scheitelpunkts der Parabel.