Genauigkeitsproblem

[Philip K. (plip)]

Hallo,

ich stehe im Moment bei meiner Diplomarbeit ein wenig auf dem 
Schlauch...

Ich habe ein PCM-Eingangssignal mit 16 Bit Integer Auflösung. Dieses 
durchläuft eine kurze FFT und ein paar FIR-Filter. Dabei werden die 
Werte natürlich zwangsläufig größer und es kommt möglicherweise zu einem 
Overflow.

Ich habe also prinzipiell folgende möglichkeiten(ohne Anspruch auf 
Vollständigkeit):

a) Ich arbeite mit Floating-Point
b) Ich skaliere die Werte runter
c) Ich vergrössere bei jeder Berechnung die Auflösung entsprechend der 
Operation

a) habe ich eigentlich zu Beginn der Arbeit ausgeschlossen, weil das 
Sytem mit möglichst wenig Ressourcen auskommen soll. Außerdem müsste ich 
mein halbes Design wieder umschmeißen...
b) Geht natürlich zu Lasten der Genauigkeit. Allerdings bin ich nicht 
auf eine durchgehend korrekte Verarbeitung angewiesen, soll heißen ich 
könnte zur Laufzeit auf Overflows reagieren und meine Skalierung 
dynamisch machen.
c) Wäre insofern reizvoll, weil ich halt keinen Genauigkeitsverlust 
habe. Und die Wortbreite ist ja auf einem FPGA quasi beliebig. 
Allerdings stoße ich dabei auf das Problem, daß ich aus komplexen 24-Bit 
Werten den Betrag berechnen muss, das bedeutet 24 Bit Multiplikation und 
49 Bit Wurzel ziehen. Ersteres ist ja noch relativ unproblematisch aber 
für zweiteres habe ich keine fertige Lösung gefunden, müsste also selbst 
per Newton oder so eine Funktion entwerfen. Die dürfte dann aber auch 
relativ Ressourcenfressend sein.

Wie würdet Ihr vorgehen?

[sechsdreizwei (Gast)]

Wenn 16 bit nicht genuegen, vielleicht 24bit, 32 bit ?

[Philip K. (plip)]

Da hab ich dann halt das Problem der Betragsbildung. 24 Bit zum Quadrat 
macht 48. 48 Bit + 48 Bit macht 49 und daraus muss ich dann die Wurzel 
ziehen. Gefunden habe ich nur eine Wurzelfunktion bis 32 Bit. Damit 
hätte ich schon bei 16 Bit Auflösung ein Problem.

[Falk B. (falk)]

@ Philip Kirchhoff (plip)

>Ich habe ein PCM-Eingangssignal mit 16 Bit Integer Auflösung. Dieses
>durchläuft eine kurze FFT und ein paar FIR-Filter. Dabei werden die
>Werte natürlich zwangsläufig größer und es kommt möglicherweise zu einem
>Overflow.

Worauf? FPGA oder DSP oder uc?

>a) Ich arbeite mit Floating-Point

Würde ich nicht machen.

>b) Ich skaliere die Werte runter

Unschön.

>c) Ich vergrössere bei jeder Berechnung die Auflösung entsprechend der
>Operation

Im FPGA meist günstig machbar.

>könnte zur Laufzeit auf Overflows reagieren und meine Skalierung
>dynamisch machen.

UHHH, wenn das mal nicht daneben geht. Klingt heiss und aufwändig.

>49 Bit Wurzel ziehen. Ersteres ist ja noch relativ unproblematisch aber
>für zweiteres habe ich keine fertige Lösung gefunden, müsste also selbst
>per Newton oder so eine Funktion entwerfen. Die dürfte dann aber auch
>relativ Ressourcenfressend sein.

Warum? Du hast IMMER die Möglichkeit, zwischn einer rein seriellen, 
kleinen, langsamen und voll parallelen, grossen, schnellen Lösung zu 
wählen. Wenn dein Datendurchsatz nicht sehr hoch ist würde ich die 
serielle Lösung wählen, die dann auch sehr ähnlich zu uC Lösungen 
aussieht.

>ziehen. Gefunden habe ich nur eine Wurzelfunktion bis 32 Bit. Damit
>hätte ich schon bei 16 Bit Auflösung ein Problem.

Die sollte man doch relativ einfach auf 49 Bit aufbohren können, oder?

MFg
Falk

[Morin (Gast)]

> a) Ich arbeite mit Floating-Point

Floating-Point passt dann, wenn du eine hohe relative Genauigkeit 
sowohl bei sehr großen als auch sehr kleinen Werten brauchst. Wenn es um 
hohe absolute Genauigkeit geht, nimm integer bzw. Fixed-Point.

Du definierst dein Eingabesignal als "16-bit integer", also einen 
Bereich von 2^16 Werten mit gleichbleibender absoluter Genauigkeit. Da 
bringt dir Floating-Point nicht viel.

> Dieses durchläuft eine kurze FFT und ein paar FIR-Filter.

Du müsstest doch eigentlich eine Vergößerung deines Wertebreichs 
berechnen können, also den größten Abtastwert, der in einem gefilterten 
Signal vorkommen kann. Wenn du das machst, bekommst du eine Antwort für 
die Breite das Ausgabeworts. Das sagt zwar noch nichts über 
Zwischenergebnisse aus, aber immerhin.

Dasselbe dann für den größten Koeffizienten nach der FFT.

[Falk B. (falk)]

@ Morin (Gast)

>Floating-Point passt dann, wenn du eine hohe relative Genauigkeit
>sowohl bei sehr großen als auch sehr kleinen Werten brauchst.

Schon, aber Floating Point braucht auf nem FPGA tierisch Resourcen.

MFG
Falk

[Günter -. (guenter)]

Philip Kirchhoff wrote:
[...]
> c) Ich vergrössere bei jeder Berechnung die Auflösung entsprechend der
> Operation
[...]
> c) Wäre insofern reizvoll, weil ich halt keinen Genauigkeitsverlust
> habe. Und die Wortbreite ist ja auf einem FPGA quasi beliebig.
> Allerdings stoße ich dabei auf das Problem, daß ich aus komplexen 24-Bit
> Werten den Betrag berechnen muss, das bedeutet 24 Bit Multiplikation und
> 49 Bit Wurzel ziehen. Ersteres ist ja noch relativ unproblematisch aber
> für zweiteres habe ich keine fertige Lösung gefunden, müsste also selbst
> per Newton oder so eine Funktion entwerfen. Die dürfte dann aber auch
> relativ Ressourcenfressend sein.
>
> Wie würdet Ihr vorgehen?

Was für eine FFT nimmst du denn? Hast du die selber entwickelt?

Bei einer radix-2 FFT werden die Daten maximal um 1-Bit pro Rank (Spalte 
mit Butterfly Operationen) wachsen, vorausgesetzt der Twiddlefaktor 
Multiplizierer skaliert das Ergebnis wieder auf die Eingangsbitweite 
zurück. Für N=256 z.B. werden die Ausgangsdaten 8-Bit größer sein als 
die Eingangsdaten.

In Abhängigkeit der Eingangsdaten ist es nicht immer nötig das volle 
Wachstum mitzumachen und ein Skalieren in selektierten Ranks ist 
möglich. Um das heraus zu finden ist es aber nötig die FFT mit genügend 
Testdaten zu simulieren die dem erwarteten Eingangssignal entsprechen.

[Philip K. (plip)]

Die FFT kommt vom Xilinx Core Generator. Aber auch die kann man manuell 
skalieren. Die Eingangswerte haben eine beliebige Größe.

Mein Problem ist ja auch nicht, dass ich nicht weiß, wie groß die Werte 
werden können, sondern daß sie eine so große Wortbreite haben, dass ich 
für Standardoperationen keine fertigen Lösungen finde. Und jedes kleine 
Teil selbst zu bauen ist halt Zeitaufwändig. Sollte aber schon noch drin 
sein...

@Falk
Was die Wurzel angeht. Leider ist der Code nicht dokumentiert und ich 
versteh nicht wirklich was er macht. Außerdem ist es eine parallele 
Lösung. Ich denke ich werd dann doch was selbst bauen. Nach allem was 
ich bis jetzt gelesen habe bietet sich das Heron-Verfahren an. 
Allerdings benötige ich dann auch einen großen Dividierer aber um den 
komm ich wohl eh nicht herum.

[Morin (Gast)]

> Schon, aber Floating Point braucht auf nem FPGA tierisch Resourcen.

Klar, aber wenn du entsprechende Anforderungen an die Genauigkeit hast 
dann hast du sie eben, Ressourcenverbrauch hin oder her.

Einer meiner Profs hat von einem Projekt berichtet, in dem ein 
MP3-Player gebaut wurde. Da wurde auch des Ressourcenverbrauchs wegen 
Fixed-Point statt Floating-Point verwendet. Ergebnis: der Klang war 
miserabel. (Details kann ich nicht erzählen da ich an dem Projekt nicht 
beteiligt war).

Gerade deshalb ist es wichtig zu wissen, ob man solche Anforderungen 
an die Genauigkeit hat oder nicht.

[Falk B. (falk)]

@ Morin (Gast)

>Klar, aber wenn du entsprechende Anforderungen an die Genauigkeit hast
>dann hast du sie eben, Ressourcenverbrauch hin oder her.

Ganz so einfach ist es nicht. Man kan mit Festkomman noch EINIGES 
machen, z.B. dynamische Anpassung der Datenbreite. geht dann in Richtign 
Flisskomma, ist es aber noch nicht. Kostet wenig mehr als Festkomma, 
kann aber fast soviel wie Fliesskomma.

MFG
Falk

[Andreas S. (andreas) (Admin)]

Wegen Fixed Point muss kein MP3-Decoder schlecht klingen, wie u.a. MAD 
beweist. Allerdings muss man ein bisschen mehr nachdenken, einfach 
runterschreiben wie in Floating Point ist nicht drin.

[Philip K. (plip)]

So, das Wurzelziehen funktioniert. Habe mich jetzt für die Fixed Point 
Variante entschieden. Dividierer und Iteration laufen komplett seriell. 
Benötigt zwar, je nach Eingangswert, bis zu 800 Zyklen bei 49 Bit, aber 
die Zeit hab ich...
Zwar könnte ich noch was drehen, wenn ich den Startwert der Iteration an 
den Eingangswert anpasse aber dafür fehlt mir die Zeit. Und optimieren 
kann man ja am Ende immer noch.

[Falk B. (falk)]

Schön das zu hören.

MfG
Falk

[my name is nobody (Gast)]

Hallo,
habe einige FFT-Routinen für Microcontroller geschrieben u.A für 8051 in 
Assembler. Habe da den Trick gefunden, nach jedem zweiten Rank (Spalte
mit Butterfly Operationen) einmal das ganze Datenfeld um ein Bit nach 
rechts zu schieben (Aritmetic shif right). So habe ich meine daten vor 
dem Überlaufen bewahrt.
Das wird bei deinem fertigen FFT wohl nicht möglich sein.

Zweiten kannst du die 48 Bit Zahlen vor der Wurzelberechnung normieren:
D.H. du schiebst soweit nach links bis das erste Bit eine "1" ist.
Anzahl der Schiebungen auf 32 limitieren.
Du merkst dir die Anzahl der Schiebevorgänge.
Von deinem Datenwort nimmst du jetzt die obersten 16 Bit und ziehst die 
Wurzel. Ergebnis 8 Bit.
Diese 8 Bit schiebst du um die halbe Anzahl der gemerkten 
Schiebevorgänge nach links.
Sind jetzt 24 Bit Wortbreite.
Falls die Anzahl der gemerkten Schiebevorgänge eine ungerade Zahl ist,
noch mit der Konstante Wurzel 2 multiplizieren.

Alles klar oder versteht das wieder keine Sau?

[Philip K. (plip)]

my name is nobody wrote:
> Hallo,
> habe einige FFT-Routinen für Microcontroller geschrieben u.A für 8051 in
> Assembler. Habe da den Trick gefunden, nach jedem zweiten Rank (Spalte
> mit Butterfly Operationen) einmal das ganze Datenfeld um ein Bit nach
> rechts zu schieben (Aritmetic shif right). So habe ich meine daten vor
> dem Überlaufen bewahrt.
> Das wird bei deinem fertigen FFT wohl nicht möglich sein.

Doch, ist möglich und wird auch gemacht. Ein fertiger FFT-Core, der 
Überlauf erzeugt macht kaum Sinn. :-)

> Zweiten kannst du die 48 Bit Zahlen vor der Wurzelberechnung normieren:
> D.H. du schiebst soweit nach links bis das erste Bit eine "1" ist.
> Anzahl der Schiebungen auf 32 limitieren.
> Du merkst dir die Anzahl der Schiebevorgänge.
> Von deinem Datenwort nimmst du jetzt die obersten 16 Bit und ziehst die
> Wurzel. Ergebnis 8 Bit.
> Diese 8 Bit schiebst du um die halbe Anzahl der gemerkten
> Schiebevorgänge nach links.
> Sind jetzt 24 Bit Wortbreite.

Ok, auf die Idee, mir die Shifts zu merken bin ich nicht gekommen. 
Trotzdem bedeutet das, je nach Wert, Datenverlust.

> Falls die Anzahl der gemerkten Schiebevorgänge eine ungerade Zahl ist,
> noch mit der Konstante Wurzel 2 multiplizieren.

Was in Integer einer 1 entspricht.

> Alles klar oder versteht das wieder keine Sau?

Doch, geht schon...;-)

[my name is nobody (Gast)]

Habe mir dazu noch folgendes überlegt:
Die Anzahl der Schiebevorgänge nach links bis die erste "1" kommt sollte 
eine gerade Zahl sein (0,2,4,6...) Damit erspart man sich die dämliche 
Multiplikation mit Wurzel 2.

Im "Schaltplan" würde ich je zwei aufeinander folgende Bits des 
Datenworts miteinander oder verknüpfen. Diese "veroderten" Bits dann auf 
einen Prioritätsdecoder.  In TTL währe das sowas wie ein 74xy147 oder 
74xy148.
Dieser steuert dann den "Barrel-Shifter"

Der Genauigkeitsverlust ist IMHO zu vernachlässigen.

[Philip K. (plip)]

Was ist ein Prioritätsdecoder?

[Philip K. (plip)]

my name is nobody wrote:

> Der Genauigkeitsverlust ist IMHO zu vernachlässigen.

Ist das eine Schätzung, Erfahrung oder hast Du es nachgerechnet?

[Walter (Gast)]

falls du mit deinem Wurzelziehen noch unzufrieden bist:
schau die mal den Cordic Algorithmus an, einer der schnellsten 
Algorithmen

[Philip K. (plip)]

Hört sich interessant an! Mich interessiert weniger die Geschwindigkeit 
als der Ressourcenverbrauch. Und der Cordic kommt angeblich nur mit 
Schieben und addieren aus. Meine Lösung dagegen benoetigt einen großen 
Dividierer.
Wenn ich am Ende noch Zeit habe, werd ich da mal dran drehen...

[Stefan Hanke (Gast)]

Im aktuellen 9.2er ISE haben die den Core Generator aber ordentlich 
beschnitten. Ich weiß, dass vorherige Versionen auch Dividierer und 
Wurzelzieher als Operation angeboten wurden, die mittels CORDIC 
implementiert wurden. Man konnte zwischen voll-parallel und serieller 
Implementierung wählen, außerdem noch "Coarse Rotation" etc.

Falls du mit dem 9.2er arbeitest, solltest du mal eine frühere Version 
(9.1er oder ne 8.2er) installieren, ein Projekt erstellen und die Cores 
anschauen, die dir der Core Generator dort anbietet. Wenn du dann die 
*.xco-Dateien rüberkopierst, klappt es hoffentlich -- probiers einfach 
mal aus.

[my name is nobody (Gast)]

>Was ist ein Prioritätsdecoder?

Eine Schaltung die eine bestimmte Anzahl von Leitungen überwacht und 
wenn eine davon auf "1" liegt wird am Ausgang binär die nummer der 
Leitung ausgegeben. Wenn mehrere der Eingangsleitungen auf "1" liegen, 
wird die mit der höchsten Priorität ausgegeben. Deshalb einfach mal 
Datenblatt vom 74HC147 oder 74HC148 anschauen.

>> Der Genauigkeitsverlust ist IMHO zu vernachlässigen.

>Ist das eine Schätzung, Erfahrung oder hast Du es nachgerechnet?

Ist Erfahrung, aber vertraue mir nicht, teste es mal selber aus, z.B. 
mit einem kleinen Basic-Wegwerfprogramm.

Währe schön, wenn du deinen Cordig-Wurzel-Algorithmus hier mal posten 
könntest. Meine Assembler-Microcontroller-Wurzel-Routine braucht einen 
Multiplizierer und probiert jedes Bit dann aus. Dauert also schon ein 
paar Zyklen...

[my name is nobody (Gast)]

Hier ein Beispiel für einen Prioritätsdecoder,
programmiert mit ORCAD PLD V1.2

|  LOW:      SW1, SW2, SW3, SW4, SW5, SW6, SW7 | Eingänge
|  HIGH:     X0, X1, X2 | Ausgänge
|
|| ---------- EQUATES ----------
|  GL7      =    SW7                     | SW7 gedrueckt
|  GL6      =    SW6  & SW7'             | SW6 & SW7 nicht gedrueckt
|  GL5      =    SW5  & SW[7~6] == 0     | SW5 & SW6-7 nicht gedrueckt
|  GL4      =    SW4  & SW[7~5] == 0     | SW4 & SW5-7 nicht gedrueckt
|  GL3      =    SW3  & SW[7~4] == 0     | SW3 & SW4-7 nicht gedrueckt
|  GL2      =    SW2  & SW[7~3] == 0     | SW2 & SW3-7 nicht gedrueckt
|  GL1      =    SW1  & SW[7~2] == 0     | SW1 & SW2-7 nicht gedrueckt
|
|  X0   =    GL1 # GL3 # GL5 # GL7
|
|  X1   =    GL2 # GL3 # GL6 # GL7
|
|  X2   =    GL4 # GL5 # GL6 # GL7

Und nun das Datenblatt vom 74HC147

[Philip K. (plip)]

Ohje, jetzt muss ich meine unausgereiften Algorithmen auch noch 
öffentlich machen...:-)

Grundlage findest Du hier: http://de.wikipedia.org/wiki/Heron-Verfahren

entity sqrtNBit_ser is
  generic(
      N: integer := 49
  );
  port(
    input    : in std_logic_vector(N-1 downto 0);
    sqrtout  : out std_logic_vector((N+1)/2-1 downto 0);
    clk    : in std_logic;
    start    : in std_logic;
    done    : out std_logic
    
  );
end sqrtNBit_ser;

architecture bhv of sqrtNBit_ser is
  signal temp      : std_logic_vector(N-1 downto 0);
  signal quotient  : std_logic_vector(N-1 downto 0) := (others => '0');
  signal sum      : std_logic_vector(N-1 downto 0) := (others => '0');
  signal start_div  : std_logic := '0';
  signal done_div  : std_logic := '0';
begin
  Divide: dividerNBit_ser
    generic map(N)
    port map(
      divident  => input,
      divisor  => temp,
      quotient => quotient,
      clk    => clk,
      start    => start_div,
      done    => done_div
    );
    
  sum <= temp + quotient;
  sqrtout <= temp((N+1)/2-1 downto 0);
    
  HeronSM: process(clk)
    type state is(idle, dividing, comparing);
    variable pstate      : state:= idle;
    variable switch_state  : boolean := false;
    variable temp_alt      : std_logic_vector(N-1 downto 0);
    variable c          : integer := 0;
    variable waitfordiv    : boolean := false;
    variable again        : boolean := false;
    variable waitcycle    : boolean := false;
  begin
    if(rising_edge(clk))then
    
      if(switch_state)then
        switch_state := false;
        waitfordiv := false;
        waitcycle := false;
        if(pstate = idle)then
          pstate := dividing;
        elsif(pstate = dividing)then
          pstate := comparing;
        elsif(pstate = comparing)then
          if(again)then
            pstate := dividing;
            again := false;
          else
            pstate := idle;
          end if;
        end if;
      end if;
      
      if(pstate = idle)then
        if(start = '1')then
          done <= '0';
          switch_state := true;
          temp <= (others => '0');
          temp_alt := (others => '0');
          temp(N/4+2) <= '1'; --Startwert
        end if;
      end if;
      
      if(pstate = dividing)then
        if(not waitfordiv)then
          --divisor <= temp;
          start_div <= '1';
          waitfordiv := true;
          waitcycle := true;
        elsif(waitcycle)then
          waitcycle := false;
        else
          start_div <= '0';
          if(done_div = '1')then
            --temp(N-2 downto 0) <= sum(N-1 downto 1);
            switch_state := true;
            --itstep <= '1';
          end if;
        end if;
      end if;
      
      if(pstate = comparing)then
      
        if(temp(N-2 downto 0) = sum(N-1 downto 1))then 
          -- Ergebnis ist vollst. konvergiert
          done <= '1';
        else
          if(temp_alt(N-2 downto 0) = sum(N-1 downto 1))then
            -- Ergebnis alterniert --> vorigen Wert verwenden
            temp(N-2 downto 0) <= temp_alt(N-2 downto 0); 
            done <= '1';
          else
            -- Weitermachen
            again := true;
            temp(N-2 downto 0) <= sum(N-1 downto 1);
          end if;
        end if;
        switch_state := true;
        temp_alt := temp;
        
      end if;
    end if;
  end process;
end bhv;

[Philip K. (plip)]

Hab jetzt noch ein Lookuptable für den Startwert vorgeschaltet und die 
Iterationszahl damit um den Faktor 3-4 reduziert.