FIR Filter um Sinus zu detektieren

[Techniker (Gast)]

Hallo,

ich möchte mit einem FIR Filter einen Sinus detektieren. Der Sinus liegt 
in einem definierten Zeitbereich. Könnte man hierfür einen FIR Bandpass 
einsetzen ?

Für weitere Ratschläge bin ich sehr dankbar!

[Goerzel (Gast)]

Für die einfache detektion von Frequenzen (z.B. DTMF) wird meist der 
Goerzelalgorithmus verwendet.
Siehe z.B. hier:
http://netwerkt.wordpress.com/2011/08/25/goertzel-filter/

[Techniker (Gast)]

Guten Morgen,

besten Dank für den Tipp. Gibt es in Matlab für den Goerzelalgorithmus 
eine Funktion ?

[Techniker (Gast)]

Der Goertzelalgorithmus ist ja im prinzip ein IIR Filter. Ich möchte 
aber einen FIR-Filter hierfür verwenden. Würde es eigentlich Sinn machen 
einen Bandpass FIR-Filter zu verwenden ?

[Techniker (Gast)]

In Matlab gibt es eine Funktion mit dem Namen goertzel. Leider besitze 
ich diese Funktion nicht. Nun habe ich mir gedacht ich programmiere den 
Goertzlalgorithmus mit der filter Funktion nach.

Siehe Übertragungsfunktion:
http://www.mathworks.de/help/toolbox/signal/ref/goertzel.html

%N = 256
for i=0:1:N-1
  f = filter([1 -exp(-i*2*pi*i/N)],[1 -2*cos(2*pi*i/N) 1],1:N);
end;

Müsste diese Implementierung so stimmen ?
Wie füge ich nun das Sinussignal (100kHz mit N = 256 Abtastwerte) hinzu 
?

[H.B. (Gast)]

Sehe ich nicht, wie das gehen soll, was macht die Filterfunktion?

Musst Du nicht einfach nur für jeden Wert die jeweils beiden 
zurückliegenden Werte aufgreifen und mit dem COS-Faktor multiplizieren, 
wie es in WIki steht? Ich meine das englische Wiki!

[Strubi (Gast)]

Tip: Etwas generischer geht's mit der Signal Processing toolbox so:

fgoer = tf([0.5], [ 1, -(r*kcos1), r*r ]);
filtered = filter(tf.num, tf.den, data);

tf plus entsprechende Argumente baut die Transfer-Funktion für den 
Goertzel.
r ist typischerweise fast 1 (0.9999...), kcos1 der cos(2*pi*k/N)-Term 
für die gewünschte Frequenz.

[Strubi (Gast)]

Da hatte ich zu schnell geklickt: In der zweiten Zeile sollte tf = 
fgoer.

Sorry.

[Techniker (Gast)]

Hallo Strubi,

danke für deine Hilfe. Ich nicht nachvollziehen wie du auf diese Zeile 
kommst:
>>fgoer = tf([0.5], [ 1, -(r*kcos1), r*r ]);

Mir fehlt eine genaue Herleitung des Goerzel Algorithmus.

[Techniker (Gast)]

Ich habe den Filter so implementiert wie auf der Internetseite 
http://www.google.de/url?sa=t&rct=j&q=sa2011-lab3-fourier&source=web&cd=1&ved=0CCcQFjAA&url=http%3A%2F%2Fwww.emt.tugraz.at%2Fsystem%2Ffiles%2FSA2011-Lab3-Fourier.pdf&ei=8y6RT8aELo2GhQfkhb2fBA&usg=AFQjCNEok5d4o0Efmwge58wWGKVdjsM_fw 
beschrieben (Seite 9).

f=filter([1],[1 -exp(-j*2*pi*k/N)],y)

y sind meine Eingangsdaten
N ist Konstant
und k ist der Laufindex von k=0... N-1

Nach der Berechnung mache ich dann eine Betragsbildung mit dem Befehl 
abs.
Stimmt meine Implemntierung in Matlab ?

[Goerzel (Gast)]

Techniker schrieb:
> Stimmt meine Implemntierung in Matlab ?

Wie wäre es wenn Du Dir einfach ein Testsignal erzeugst, z.B. so was:

f=50;            %Signalfrequenz
fs=100*f;        %Abtastfrequenz
t=0:1/fs:1/f;    %Zeitvektor erstellen

y= sin(2*pi*f*t)+sin(2*pi*f*t)+sin(4*pi*f*t);

Dann lässt Du den Filter drüber laufen.
Anschließend ne FFT und Du weist es.

Grüße

[Techniker (Gast)]

Irgendwie funktioniert das ganze nicht so richtig.
Was mache ich da noch falsch ?

f=50;            %Signalfrequenz
fs=100*f;        %Abtastfrequenz
t=0:1/fs:1/f;    %Zeitvektor erstellen

y= sin(2*pi*f*t);
max_y = max(abs(y))*1.1; 
subplot(2,1,1);
fig = figure(1); 
plot(t,y);
axis([0 t(end) -max_y max_y]);
grid on;

N = 128;      % Abtastwerte
k = 0:1:N-1;
g=filter([1],[1 -exp(-j*2*pi*k/N)],y);

H = fft(g,N); 
amplH = abs(H);
amplitudengang = fftshift(amplH/N);
subplot(2,1,2);
stem(amplitudengang);
grid on;

[Mathegenie (Gast)]

Was funktioniert nicht richtig? Siehst du Deine Frequenz nicht, oder 
siehst Du gar nichts?

[J. S. (engineer)]

Mathegenie schrieb:
> Was funktioniert nicht richtig? Siehst du Deine Frequenz nicht, oder
> siehst Du gar nichts?

Kann es sein, dass Du Deine Frequenz deshalb nicht findest, weil Du zu 
grob abtastest?

Die Schrittweite der Testfrequenzen mit denen Du rechnest, muss in einem 
sinnvollen Verhältnis zur Samplezahl stehen, weil mit zunehmender 
Samplezahl die DFT-Rechnung immer selektiver wird und Dir so eine 
Frequenz durchrutscht. Ich habe im Anhang ein altes Beispiel hinkopiert:

Die Schrittweite von 8 führt dazu, dass bei der blauen Kurve mit den 
Punkten die Frequenz einigermassen getroffen wird und sich ein peak 
bildet, bei 10 Schritten liegt man haarscharf daneben und der peak geht 
im Rauschen unter. Bei der gröber auflösenden, türkisen Kurve wurde die 
Berechung schon nach 2048 von 4096 samples abgebrochen.

Umgekehrt führt die Überlegung dahin, dass Du Deine Sanplezahl nicht zu 
gross wählen darfst, wenn Du Deine Frequenz nicht genau triffst, oder Du 
bekommst eben das Ergebnis, daß das Matching gering ist.


Techniker schrieb:
> Würde es eigentlich Sinn machen
> einen Bandpass FIR-Filter zu verwenden ?
Der Bandpass ist nicht anderes als eine Summe von Frequenzen als Maske, 
die gleichzeitig auf den Datenstrom losgelassen werden. Mathematisch ist 
der Bandpass ein Hochpass, von dem ein TP abgzogen wurde. In den 
Koeffizienten zeigt sich das durch ein Integral der oberen Grenzfrequenz 
hin zu 0, von der der Koeefizientensatz unterhalb der GF weggelassen 
wurde.

Wenn Du dir einander annäherst, hast Du am Ende eine Frequenz wie oben.

Wenn Du Deine Frequenz nicht genau kennst, ist der BP gfs besser.

[Techniker (Gast)]

Hallo Jürgen,

danke für deinen Beitrag. Du meinst ich taste zu grob ab. Auf einem 
System kann habe ich nun mal pro Sinusperiode 128 Abtastwerte. Das 
müsste doch reichen. Ich habe die Vermutung das irgendwo noch ein Fehler 
drin ist. Bin leider noch nicht weitergekommen.

[Techniker (Gast)]

Mit dem Goertzelalgorithmus möchte ich erkennen, ob der Sinus zum 
Beispiel mit 100Khz vorhanden ist oder nicht. Wenn Impulsstörungen oder 
auch Rauschen additiv mit dem Sinus überlagert sind, soll trotzdem der 
Sinus von 100Khz detektiert werden.

[tsag (Gast)]

Techniker schrieb:
> vorhanden ist oder nicht
ist eine Frage der Definition. Du wirst beim Rauschen immer einen mehr 
oder weniger grossen Wert für die Korrelation erhalten. Du brauchst dann 
einen Grenzwert zum Entscheiden.

[J. S. (engineer)]

Techniker schrieb:
> danke für deinen Beitrag. Du meinst ich taste zu grob ab. Auf einem
> System kann habe ich nun mal pro Sinusperiode 128 Abtastwerte.

Ich denke, wenn dann diesbezüglich "zu fein". Je mehr Punkte je Sinus, 
desto selektiver wird diese Maske, mit der Du den Datenstrom ansiehst. 
Wenn diese Punktezahl steigt, muss auch die "Maskenanzahl" steigen.

In der Grafik habe ich das veranschaulicht:

* Blau: Sehr genau abgetastet, man erkennt den Verlauf des Signals
* Rot: Mit einem 7-tel der Aufösung abtastet: Das Maximum rutscht durch
* Grün: Mit einem 7-tel der Auflösung bei 7-facher Breite
* Hellgrün genau abgetastet mit 7facher Breite

Der grüne ist also grob/breit genug, um das Signal noch zu erfassen und 
fein genug, um den Verlauf noch erahnen zu können.

Blau wäre etwa die obere Grenze. Grün etwa die untere Grenze.
.