FIR in Octave

[U. B. (ub007)]

Hallo !

Es gab schon mal einen ähnlichen Thread hier aber genau auf die Frage wo 
ich jetzt auch habe wurde leider keine Antwort gegeben.
Ich möchte in Octave einen FIR berechnen aber nicht mit den Funktionen 
die Octave zur Verfügung stellt sondern eben "mit der Hand am Arm".
Ich habe eine Sweep-WAV generiert die ich in Octave einlesen kann. Die 
Frequenz läuft von 50Hz bis 18kHz mit 22050 Samples bei einer Samplerate 
von 44100Hz. Mein LP-Filter soll z.B. bei 2.5KHz dicht machen. Wie 
erwartet zeigt die FFT-Analyse (siehe Bild) wie der Sweep im 
Frequenz-Bereich aussieht. Lass ich aber den Filter mit 20 Tabs drüber 
laufen ergibt sich nur in der Amplitude eine Änderung. Auch mit 100 Tabs 
oder ungeraden oder geraden Anzahl von Tabs sieht das immer gleich aus.
OK, der Code in Octave ist nicht effizient, aber wie gesagt, es sollte 
ja erst mal nur zum testen sein. Die Koeffizienten wurden über Tone FIR 
Designer errechnet und der Code habe ich von
http://www.fh-schmalkalden.de/schmalkaldenmedia/Realisierung_Digitaler_Filter_in_C-p-419.pdf
An was kann es liegen das man ab 2.5kHz keinen Abfall der Frequenz sieht 
?

Gruß Uli

Hier mal der Octave-Code:

% Octave  

%// FilterType = LoPass, fcut1 = 2.5kHz 
%// fsample = 44100Hz, Window = Hamming, 
nc = 20;
b=[
  -0.011433, 0.007947, 0.060518, 0.184715, 0.406822, 0.726661,
  1.111126, 1.498814, 1.815087, 1.992966, 2.003541, 1.815087,
  1.498814, 1.111126, 0.726661, 0.406822, 0.184715, 0.060518,
  0.007947, -0.011433,
];

% ---------------------------------------- 
function y=G(ys,nc,b,in_buffer) 

   % alle Werte nach rechts shiften
  for i=nc:-1:2          
      in_buffer(i) = in_buffer(i-1); 
   end 

   % Schreibe neuen Wert ys in Buffer 
   in_buffer(1) = ys; 
   
   % Berechne neuen Ausgangswert 
   y = 0; 
   for i=1:+1:nc 
      y += ( b(i)*in_buffer(i) ); 
   end 
endfunction 
% ---------------------------------------- 
clf
t=1; y=0; ys1=0; samplerate=44100;
ys = wavread ("ttt.wav");
samples=length(ys);

% in_buffer loeschen
for i=1:+1:nc
   in_buffer(i) = 0; 
end

for t=1:+1:samples
   ym(t)=G(ys(t),nc,b,in_buffer);
end 

ym1=abs(fft(ym));
yo1=abs(fft(ys));

t=1:+1:samples;

% FFT-Plot
subplot(4,1,1)
plot(t*(samplerate/samples),yo1);grid on;
axis([0 20000])
title("FFT ohne Filter")

subplot(4,1,2)
plot(t*(samplerate/samples),ym1);grid on;
axis([0 20000])
title("FFT nach Filter")

%Signal-Plot
subplot(4,1,3)
plot(t,ys(t));grid on;
title("Signal ohne Filter")

subplot(4,1,4)
plot(t,ym(t));grid on;
title("Signal mit Filter")

[Markus B. (russenbaer)]

Servus Uli,

Ich hab Dir mit Octave (freqz Funktion) den Frequenzgang Deines Filters 
gerechnet.
Willst Du wirklich so ein Filter?

BTW:
Du gibst den State (Gedächtnis) Deines Filters nicht zurück, MATLAB und 
ich nehme an auch OCTAVE übergeben Werte durch kopieren, nicht durch 
Referenz. Im (hab ich mir jetzt nicht angeschaut) Beispiel wo Du den 
C-Code nachgetippt hast wird sicher ein Pointer auf deinen in_buffer 
übergeben.
Das gibts in MATLAB/Octave nicht!

Du musst als zweite Variable den State zurückgeben!

function [y, in_buffer] = G(ys,nc,b,in_buffer)

lg
Markus

PS:
Effizientes umkopieren (was Du in einer Schleife erledigst)

b(2:end) = b(1:end-1];

Inprodukt zweier Vektoren (Berechnung Deines Ausgangswertes y, Ersatz 
für die zweite Schleife in Deiner Funktion)

y = b(:)' * in_buffer(:);

Erzeugen eines Null-Initialisierten Buffers:

in_buffer = zeros(numel(b),1);

[U. B. (ub007)]

Hallo Markus !

Erstmal danke für die Antwort.

Markus B. schrieb:
> Servus Uli,
>
> Ich hab Dir mit Octave (freqz Funktion) den Frequenzgang Deines Filters
> gerechnet.
> Willst Du wirklich so ein Filter?
>
Nein Markus natürlich nicht, aber das ist wie beschrieben eben nur zum 
Test.
>
> BTW:
> Du gibst den State (Gedächtnis) Deines Filters nicht zurück, MATLAB und
> ich nehme an auch OCTAVE übergeben Werte durch kopieren, nicht durch
> Referenz. Im (hab ich mir jetzt nicht angeschaut) Beispiel wo Du den
> C-Code nachgetippt hast wird sicher ein Pointer auf deinen in_buffer
> übergeben.
> Das gibts in MATLAB/Octave nicht!
>
Ebent ;-) Das gibt es nicht. Deshalb hab ich mal abgeprüft ob Variablen 
wie in_buffer oder eben die Koeffizienten b(i) in der Funktion zur 
Verfügung stehen. Ich übergebe sie ja beim Aufruf mit z.B. G(...b,,,)
Und sie stehen auch drin !
>
> Du musst als zweite Variable den State zurückgeben!function [y,
> in_buffer] = G(ys,nc,b,in_buffer)
> lg
> Markus
>
Klar, man hätte korrektehalber y und in_buffer mit einem "return" 
zurückgeben können, aber nachdem die Werte außerhalb und innerhalb der 
Funktion(Scope) identisch sind(hab ich mir durch printf anzeigen 
lassen), hab ich mir das an der Stelle erspart.
Ich weiß das klingt wie hingeschludert, aber ich befürchte selbst wenn 
alles 1000% korrekt programmiert ist, dass das Ergebnis nicht besser 
aussieht als wie oben im Bild ! In dem anderen Thread sagte der 
Thread-Eröffner man muss mit den Koeffizienten einbischen 
herum-experimentieren dann bekommt man den Filter tatsächlich zum 
laufen, aber genau das möchte ich eben nicht.

Gruß Uli

[Detlef _. (detlef_a)]

Hallo Uli,

das rechtsshiften der Werte kommt nicht zurück. Machs' so wie Markus das 
beschrieben hat, Deine Variante funktioniert nicht.

Cheers
Detlef

[U. B. (ub007)]

Detlef _a schrieb:
> Hallo Uli,
>
> das rechtsshiften der Werte kommt nicht zurück. Machs' so wie Markus das
> beschrieben hat, Deine Variante funktioniert nicht.
>
> Cheers
> Detlef

Ups....da hab ich etwas nicht in Octave nicht gewußt....Shit happens...
denn als ich das von Markus einbaute...sieht das ganze jetzt so 
aus...siehe Bild...und genau das wollte ich !
Dank euch beiden :-)

Grüßle Uli