Octave/Matlab vs. GCC - sinus ist nicht gleich

[Daniel (Gast)]

Hallo,
ich habe das Problem, dass ich ein bestimmtes FFT-Fenster in
Matlab/Octave errechnen will.
Ich habe hier bereits eine Implementierung in C. Da ich
bei der Fensterung zwischen C-Code und Matlab-Code immer
noch Differenzen von 10^-12 hatte hab ich angefangen
den Fehler zu suchen.

Jetzt hab ich folgende zwei Codes-Schnipsel:
MATLAB:

t=0;
L=16000;
windowed = zeros(1,L);

while(t<L)
  wind = cos(t);
  windowed(t+1)= wind;
  t+=1;
end

C-CODE

N=16000;
int n;
float sum=0;
for (n = 0; n < N; n++)
    {
              sum = cos((float)n);
        outp[n] = sum;
    }
...
//Ausgabe:
for(i=0; i< L; i++) {
    printf("%e\n", outp[i]);
  }

So, jetzt speicher ich die GCC Werte mit ./meinprog > test.dat
und lade diese in octave rein.

Dann bilde ich die Differenz zwischen Octave-Werten (siehe ersten Code)
und GCC-Werten (siehe zweiten Code).

Im Anhang ist das Ergebnis als Plot. Wo bekommt man Infos
über die Implementierung vom Sinus oder den trigonometrischen
Funktionen in Octave und in GCC? Welche Implementierung ist jetzt
die genauere?
Kann ich von Octave auf die GCC-Sinus-Implementierung zugreifen?

Viele Grüße
 Daniel

[Andreas S. (andreas) (Admin)]

Verwende mal double statt float in deinem C-Code. Octave arbeitet 
standardmäßig mit double.

[Daniel (Gast)]

oh, da hätte ich drauf kommen können :/

Danke! Hat geholfen

[Michael L. (Gast)]

Hallo Daniel,

> oh, da hätte ich drauf kommen können :/
>
> Danke! Hat geholfen

Es gibt für Matlab eine kokstenlose Toolbox, bei der Du eine beliebige 
Genauigkeit einstellen kannst. Ich kann Dir bei Interesse den Link 
raussuchen und die Implementierung für die FFT, die ich einprogrammiert 
habe, beilegen.

Gruß,
  Michael

[Flo (Gast)]

@Michael Wäre super wenn du das mal posten könntest. Würde mich auch 
interessieren.
Gruss, Florian

[Michael L. (Gast)]

Hallo Florian,

die FFT steht unten. Vorher mußt Du die Multiple Precision Toolbox 
einbauen.
http://www.mathworks.de/matlabcentral/fileexchange/6446

Die Installation unter Windows ist nur mittelprächtig intuitiv. Du 
solltest zunächst die dll-Dateien (/toolbox/NewClasses/@mp/private) an 
eine Stelle kopieren, an der Windows bzw. Matlab sie findet und 
anschließend den Suchpfad in Matlab um den Pfad ergänzen, der die 
.m-Dateien enthält.

Folgende Dateistruktur sollte funktionieren:
Das Verzeichnis
\Matlab\toolbox\mptoolbox_1.1\
enthält verschiedene .m-Dateien (angefangen mit mp_dev_sort_bubble_a.m 
und euler.m)

in das Verzeichnis
\Matlab\toolbox\mptoolbox_1.1\@mp
kommen die restlichen .m-Dateien wie abs.m, acos.m usw

und in das Verzeichnis
\Matlab\toolbox\mptoolbox_1.1\@mp\private
kommen schließlich die zugehörigen dll- und c-Dateien.

Gruß,
  Michael



% Fouriertransformation mit beliebiger Genauigkeit
% y = fft_mp(x, prec)
%
% erlaubt Eingaben aus der multiple-precision-Toolbox.
% Die Anzahl N der Werte muß ein Vielfaches von 2 sein.
% Die Genauigkeit kann über "prec" festgelegt werden.

function ret = fft_mp(vec,prec)
  PI = 2*asin(mp(1,prec));             % pi mit hoher Genauigkeit
  n = length(vec);                     % Anzahl der Zahlen
  if (log2(n) ~= round(log2(n))),      % Fehlermeldung falls nicht 2^x
      fprintf('The number of elements must be a power of 2.\n');
      fprintf('Calculation terminated.\n');
      ret = 0;
      return;
  end;

  if n==1,                              % Rekursionsanfang
      ret = vec;
  else
      g = fft_mp(vec(1:2:n-1), prec);   % FFT-Algorithmus nach Cooley u. 
Tuckey
      u = fft_mp(vec(2:2:n),   prec);
      for k=0:(n/2-1),
          c(k+1)     = g(k+1)+u(k+1)*exp(-2*PI*i*k/n);
          c(k+n/2+1) = g(k+1)-u(k+1)*exp(-2*PI*i*k/n);
      end;
      ret = c;
end

[Thomas (Gast)]

Die Schleifen solltest du dir in Matlab & Co. tunlichst abgewöhnen.
So geht es deutlich schneller (Rechenzeit und Tipparbeit):

L = 16000;
t = 0 : L-1;
windowed = cos(t);

oder direkt:

windowed = cos(0:15999);

oder, für höhere Auflösung:

L = 16000;
Ts = 0.1;
t = 0 : Ts : L-1;
windowed = cos(t);