www.mikrocontroller.net

Forum: Digitale Signalverarbeitung / DSP Octave/Matlab vs. GCC - sinus ist nicht gleich


Autor: Daniel (Gast)
Datum:
Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
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

: Verschoben durch Admin
Autor: Andreas Schwarz (andreas) (Admin) Benutzerseite Flattr this
Datum:

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

Autor: Daniel (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
oh, da hätte ich drauf kommen können :/

Danke! Hat geholfen

Autor: Michael Lenz (hochbett)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Flo (Gast)
Datum:

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

Autor: Michael Lenz (hochbett)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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

Autor: Thomas (Gast)
Datum:

Bewertung
0 lesenswert
nicht lesenswert
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);

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.
Hinweis: der ursprüngliche Beitrag ist mehr als 6 Monate alt.
Bitte hier nur auf die ursprüngliche Frage antworten,
für neue Fragen einen neuen Beitrag erstellen.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.