---

Først og fremmest: Jeg vil ikke downloade en færdigskrevet
FFT-algoritme; jeg vil lave mig egen, fordi det er sjovest.
Først når min egen kører, vil jeg lede efter optimerede
løsninger. Jeg skriver dette, fordi jeg er træt af at få
stukket URL'er i hånden; jeg leder efter lidt grundlæggende
hjælp i matematik og programmering. Det er iøvrigt heller
ikke en FFT, jeg vil lave, men en almindelig, kedelig
integration.

I skolen skrev i et program, der kunne bestemme effekten af
frekvensindholdet. Det tog kvadratet af den diskrete Fourier-
transformation (DFT), og på den måde forsvandt alt det
imaginære, hvilket var fint nok til det, vi skulle bruge det
til (beregne forvrængning, effekt i harmoniske m.m.).

Jeg nåede aldrig helt at forstå, hvad Hanning-vinduet skulle
gøre godt for, andet end vægtningen tydeligvis gjorde mirakler,
hvis man så på signalet i gnuplot.

int
ft ( const float* x, const int N, float* X, const int K )
{
float re, im;
int n, k;


cout << "Hanning-windowing " << flush;
for ( n = 0 ; n < N ; n++ ) {
x[n] = x[n] * ( 1 - cos( n * 2*PI/(N-1) ) );
if ( n%100 == 0 ) cout << "." << flush;
}
cout << " done." << endl;


cout << "fourier-transforming " << flush;
for ( k = 0 ; k < K ; k++ ) {
re = im = X[k] = 0;

for ( n = 0 ; n < N ; n++ ) {
re += x[n] * cos( (2*PI/N)*k*n );
im += x[n] * sin( (2*PI/N)*k*n );
}

// Effektspektret: re^2 + im^2
X[k] = re*re + im*im;

if ( k % (K/10) == 0 ) cout << 100*k/K << "% " << flush;
}
cout << "100% done." << endl;


return k;
}



Men ved at tage kvadratet af talparene i sin DFT, og gøre dem
reelle, ødelægges muligheden for senere at tage en Invert
DFT, idet fortegnet ikke kan genskabet. Er det ikke rigtigt?

Nu prøvet jeg så at lave algorimen om, så den i stedet for
at gemme resultatet i en enkelt float, så gemmer den det i en
struct bestående to doubles, im og re. Og i stedet for et alm.
array, så arbejder den med en STL vector. Og vigtigst af alt:
den beregner ikke effekten.


struct complex {
double re;
double im;
};


int
fourier :: dft ( vector<complex>& x, vector<complex>& X )
{
unsigned long int n, k, N, K;
complex c;


X.clear(); // Clear any existing vector before pushing data into it.
N = x.size();
K = (N-1) / 2;

if ( verbose ) cout << "Fourier-transforming: " << flush;
for ( k = 0 ; k < K ; k++ ) {
c.re = c.im = 0;

for ( n = 0 ; n < N ; n++ ) {
c.re += x[n].re * cos( ( 2*PI / N) * k * n );
c.im += x[n].im * -sin( ( 2*PI / N) * k * n );
}

X.push_back( c );

if ( verbose ) p(k,K);
}
if ( verbose ) cout << "done." << endl;


return 0;
}


void
fourier :: hanning ( vector<complex>& x )
{
unsigned long int n, N;


if ( verbose ) cout << "Hanning-windowing: " << flush;
N = x.size();
for ( n = 0 ; n < N ; n++ ) {
x[n].re = x[n].re * ( 1 - cos( n * 2*PI/(N-1) ) );
x[n].im = x[n].im * ( 1 - cos( n * 2*PI/(N-1) ) );
if ( verbose ) p(n,N);
}
if ( verbose ) cout << "done." << endl;
}


Men et eller andet går galt. Og jeg ved ikke hvad. Jeg har
lavet to funktioner til at "konvertere" fra relle tal til
komplekse talpar og omvendt. Det er jeg nød til, da alle
beregninger skal foregå i det komplekse plan, og alt
"I/O" foregår i relle tal, f.eks. når der loades samples og
plottes.

Noget får gang og jeg kan ikke overskue det, måske nogen kan
hjælpe =).

Jeg har lagt hele herligheden ud på:
http://stephan.tisprut.dk/Skrammel/fourier/
og
http://stephan.tisprut.dk/Skrammel/fourier/fourier.tar.bz2



--
Stephan Henningsen /
/ http://tisprut.dk

---

On Wed, 11 Apr 2001 16:56:16 +0000, stephan@sanxion.tisprut.dk
(Stephan Henningsen) wrote:

>Først og fremmest: Jeg vil ikke downloade en færdigskrevet
>FFT-algoritme;

God ide, det er det man lærer mest af....

> jeg vil lave mig egen, fordi det er sjovest.
>Først når min egen kører, vil jeg lede efter optimerede
>løsninger. Jeg skriver dette, fordi jeg er træt af at få
>stukket URL'er i hånden; jeg leder efter lidt grundlæggende
>hjælp i matematik og programmering. Det er iøvrigt heller
>ikke en FFT, jeg vil lave, men en almindelig, kedelig
>integration.
>
>I skolen skrev i et program, der kunne bestemme effekten af
>frekvensindholdet. Det tog kvadratet af den diskrete Fourier-
>transformation (DFT), og på den måde forsvandt alt det
>imaginære, hvilket var fint nok til det, vi skulle bruge det
>til (beregne forvrængning, effekt i harmoniske m.m.).
>
>Jeg nåede aldrig helt at forstå, hvad Hanning-vinduet skulle
>gøre godt for, andet end vægtningen tydeligvis gjorde mirakler,
>hvis man så på signalet i gnuplot.

Vinduet "foldes" på det samplede signal. Dette gøres for at kompensere
for det faktum at du arbejder med et endeligt antal samples. Det er
ikke nødvendigt at folde et vindue på en ren samplet sinus/cosinus
kurve, hvis du blot har samplet lige præcis 1 periode af signalet. Men
hvis du derimod vil sample et signal som indeholder mange forskellige
frekvens komponenter, ja så vil du få en pludselig overgang i signalet
til den tid hvor du starter din sampling og til den tid du stopper din
sampling. Det du egentlig gør er at multiplicere det samplede signal
med et lavfrekvent signal (Hanning, Hamming, Kaiser eller lign.)

Der sker nemlig det når du sampler et signal, at du blot omsætter den
samplede "stump" til et laaaaangt repeterende signal. Jeg ved ikke om
du kan følge mig her.?

>Men ved at tage kvadratet af talparene i sin DFT, og gøre dem
>reelle, ødelægges muligheden for senere at tage en Invert
>DFT, idet fortegnet ikke kan genskabet. Er det ikke rigtigt?

Jo det er rigtig at du får effekten af den pågældende frekvens
komponent ved at kvadrere Fourier koefficenten.

Du skal først finde modulus til den pågældende frek. komponent. Det
gøres ved at finde længden af den vektor som Xre(k) + j*Xim(k)
beskriver - Pythagoras... old boy.

|X(k)| = sqrt(Xre(k)^2+Xim(k)^2)

Derefter kvadrere du så |X(k)|... for finde effekten
(Ohms lov P = U*I = U^2*R , R=1 Ohm)

Så det hele bliver faktisk...

G(k) = |X(k)|^2 = Xre(k)^2+Xim(k)^2

/Thomas

---

On Wed, 11 Apr 2001 18:52:58 +0200, Thomas Lykkeberg wrote:
>
>Der sker nemlig det når du sampler et signal, at du blot omsætter den
>samplede "stump" til et laaaaangt repeterende signal. Jeg ved ikke om
>du kan følge mig her.?

Jeg er nogenlunde med. Så for at kompenerer for, at det
ikke er et repeterende signal man sampler, men et endeligt,
hvilket giver visser fejl i udregningen, så folder man en
slags nødløsning på?


>Så det hele bliver faktisk...
>
> G(k) = |X(k)|^2 = Xre(k)^2+Xim(k)^2

Rigtigt. Og det virker fint, og det var det vores opgave
gik ud på. Men som jeg skrev, så kan man nu ikke længere
invers F-transformerer G(k), signalet kan med andre ord ikke
genskabes efter denne mishandling. Derfor vil jeg prøve at
gøre det rigtigt, om man vil. Altså ikke finde effekten,
men amplituden. Og gemme alle Fourier-koefficienter i
henholdsvis reelle og imaginære arrays.

Men det virker ikke, og jeg ved ikke hvorfor.


--
Stephan Henningsen /
/ http://tisprut.dk

---

Stephan Henningsen wrote:
> Jeg er nogenlunde med. Så for at kompenerer for, at det
> ikke er et repeterende signal man sampler, men et endeligt,
> hvilket giver visser fejl i udregningen

Hvis du tager dit signal, og repetere det så vil der være et
knæk, hvor det gentages.

Et knæk/en diskontinuitet vil give et biddrag for meget høje
frekvenser, fordi Fourier transformation er en teknik der søger
at repræsentere en diskontinuitet med vægtede kombinationer af
cos/sin funktioner med forskellige frekvens. Og det kræver en
kombination stor vægtning af hurtigt varierende funktioner,
dvs. kombinationer af høje Fourier komponenter med store amplituder
for at repræsentere en diskontinuitet.

(dette er det samme som Thomas Lykkeberg siger.)

> så folder man en slags nødløsning på?

Yeps. 1-cos(x/pi) er 0 i begge ender (0 og 1), så ganger du dit
originale signal med denne funktion så starter og slutter den
altid med 0 og er 1 stort set alle andre steder, og du kan derfor
repetere signalet uden diskontinuiteten. Du kunne lige så godt
anvende 4*x(1-x), eller (1-cos(x/pi))^n funktion der for n->oo
vil konvergere mod en funktion der er 1 over alt borset fra 0 og 1.

Resultatet er at du fjerner den støj der skyldes diskontinuiteten,
og dette bør være synligt størrelsen af amplituden af for høje Fourier
komponenter. En alternativ forklaring er at Fourier transformere
du dit oprindelige signal, så ganger du en funktion på Fourier
komponenterne der reducere de højre amplituder.

> Men som jeg skrev, så kan man nu ikke længere invers
> F-transformerer G(k), signalet kan med andre ord ikke
> genskabes efter denne mishandling. Derfor vil jeg prøve at
> gøre det rigtigt, om man vil.

Når du tager produktet så smider du fase informationen bort
for at få et reelt G(k) værdi. Uden fasen kan du ikke få det
oprindelige signal tilbage.


Btw. skal du udregne cos(n*x),sin(n*y) hurtigt så kan du bruge
additionsligningen:

cos(A+B)=cos(A)*cos(B)-sin(A)*sin(B)
sin(A+B)=sin(A)*cos(B)+cos(A)*sin(B)

med A=x og B=(n-1)x for n>=2 du kan derfor nøjes med at udregne
cos(x),sin(x) og så anvende additionsligningerne og produkter af
kendte sin og cos faktorer, dette er ideen med FFT, det smarte
i FFT er at man ikke bruger så meget hukommelse på at gemme
tabellen af sin og cos værdier.

--
* Kurosawa: http://www.designlabs.dk/husetsbio *
* Email: Carsten dot Svaneborg at risoe dot dk *
* http://www.fys.risoe.dk/fys/External/casv/ *

---

On Fri, 13 Apr 2001 18:59:56 +0200, Carsten Svaneborg wrote:
>Btw. skal du udregne cos(n*x),sin(n*y) hurtigt så kan du bruge
>additionsligningen:
>
> cos(A+B)=cos(A)*cos(B)-sin(A)*sin(B)
> sin(A+B)=sin(A)*cos(B)+cos(A)*sin(B)

Jeg har planer om at implementere en kvadratur-oscillator
som selv laver en sinus. På den måde behøver jeg heller
ikke kalde sin/cos-funktioner.

Tak for den lange forklaring. Det hjalp på forståelsen.
Det viste sig dog mit DFT-problem var en mindre og meget
dum fejl, så nu er den kørende =).

Men jeg problemer med min IDFT. Skal en invers
F-transformation ikke genskabe signalet, som det så ud, før
det blev F-transformeret? Skal man vægte det med en invers
Hanning?

Jeg har lagt mit skrammel ud. Det skulle bare være til at
pakke ud og køre 'make plot'. Så vil det (vha. gnuplot)
vise hhv. det rå signal (buffer.dec) og det samme signal
udsat for DFT og bagefter IDFT (buffer.dft). Der må være et
eller anden galt med min IDFT, måske kan du eller nogen
andre hjælpe?

http://stephan.tisprut.dk/Skrammel/fourier.tar.gz
http://stephan.tisprut.dk/Skrammel/fourier/
http://stephan.tisprut.dk/Skrammel/fourier/plot.jpg



--
Stephan Henningsen /
/ http://tisprut.dk

---

On Sun, 15 Apr 2001 22:13:41 +0000, stephan@sanxion.tisprut.dk
(Stephan Henningsen) wrote:

>Jeg har planer om at implementere en kvadratur-oscillator
>som selv laver en sinus. På den måde behøver jeg heller
>ikke kalde sin/cos-funktioner.

Jeg har før lavet en kvadratur oscillator vha. Goertzel. Den er meget
nem at implementerer, jeg ved ikke om det er nemmere end det du har
gang i.?

>Men jeg problemer med min IDFT. Skal en invers
>F-transformation ikke genskabe signalet, som det så ud, før
>det blev F-transformeret? Skal man vægte det med en invers
>Hanning?

Du har jo "filtreret" signalet ved at folde vinduet på før DFT'en, så
når du så laver din IDFT, ja så kommer du jo tilbage til signalet
efter "filtreringen"..

/Thomas

---

On Mon, 16 Apr 2001 09:38:08 +0200, Thomas Lykkeberg wrote:
>On Sun, 15 Apr 2001 22:13:41 +0000, stephan@sanxion.tisprut.dk
>(Stephan Henningsen) wrote:
>
>>Jeg har planer om at implementere en kvadratur-oscillator
>>som selv laver en sinus. På den måde behøver jeg heller
>>ikke kalde sin/cos-funktioner.
>
>Jeg har før lavet en kvadratur oscillator vha. Goertzel. Den er meget
>nem at implementerer, jeg ved ikke om det er nemmere end det du har
>gang i.?

Jeg ved faktisk ikke hvad den, jeg har lavet, hedder. Men
den ser sådan her ud:

int n, // Sampletæller.
N, // Antal samples der skal genereres
W; // Antallet af delays til beregning af w.
double phi, // Konstant beregnet på f, fs og PI.
cosphi, // Konstant, cos(phi).
sinphi, // Konstant, sin(phi).
f, // Ønskede frekvens
fs; // Samplingsfrekvens
double* w; // Mellemregninger, w[n], w[n-1] ...
double* y; // Output.
double* sinbuf; // Output.
double* cosbuf; // Output.
double x; // Input -- en delta i dette tilfælde.


N = 200;
W = 3;
f = 1000;
fs = 32000;
x = 1;

y = new double[N];
sinbuf = new double[N];
cosbuf = new double[N];

w = new double[W];
for ( n = 0 ; n < W ; n++ )
w[W-n] = 0;

phi = 2*PI*(f/fs);
cosphi = cos(phi);
sinphi = sin(phi);

for ( n = 0 ; n < N ; n++ ) {

w[W-2] = w[W-1];
w[W-1] = w[W-0];
w[W-0] = 2*cosphi * w[W-1] - w[W-2] + x;

sinbuf[n] = w[W-0] - cosphi*w[W-1];
cosbuf[n] = sinphi * w[W-1];

x = 0;
}



>>Men jeg problemer med min IDFT. Skal en invers
>>F-transformation ikke genskabe signalet, som det så ud, før
>>det blev F-transformeret? Skal man vægte det med en invers
>>Hanning?
>
>Du har jo "filtreret" signalet ved at folde vinduet på før DFT'en, så
>når du så laver din IDFT, ja så kommer du jo tilbage til signalet
>efter "filtreringen"..

Skal man så folde en invers Hanning på før at kunne
rekonstruere signalet korrekt igen? Eller er der noget
andet galt med min idft:



int
fourier :: idft ( const vector<complex>& X, vector<complex>& x )
{
unsigned long int k, K, n, N;
complex c;


x.clear(); // Clear any existing vector before pushing data into it.
K = X.size();
N = K*2;
for ( n = 0 ; n < N ; n++ ) {
c.re = c.im = 0;

for ( k = 0 ; k < N ; k++ ) {
c.re += (X[k].re + X[k].im) * cos( (2*PI/N)*k*n );
c.im += (X[k].re + X[k].im) * sin( (2*PI/N)*k*n );
}

c.re /= N;
c.im /= N;

x.push_back( c );
}


return 0;
}


complex er min egen klasse som har en double re, double im.



--
Stephan Henningsen /
/ http://tisprut.dk

---

On Mon, 16 Apr 2001 16:08:11 +0000, stephan@sanxion.tisprut.dk
(Stephan Henningsen) wrote:

>Jeg ved faktisk ikke hvad den, jeg har lavet, hedder. Men
>den ser sådan her ud:

Det er Goertzel du benytter, meget simpel kvadratur oscillator..

/Thomas

---

stephan@sanxion.tisprut.dk (Stephan Henningsen) wrote:

>Skal man så folde en invers Hanning på før at kunne
>rekonstruere signalet korrekt igen? Eller er der noget
>andet galt med min idft:
Nej det er nok ikke nogen god ide. Hvis du kun har en blok skal du
sørge for at signalet kan være i blokken uden at nå endepunkterne, så
behøver du ingen extra vinduefunktion og så kan signalet gendannes med
idft.
Hvis du derimod har delt dit signal op i flere blokke (sådan som
diverse audiofiltre virker) så skal du så vidt jeg ved lade blokkene
overlappe. Vinduefunktionen skal så være således at når du summere
blokkene bagefter giver det 1. Det er sikkert en god ide at bruge
vinduefunktionen både før og efter: vindue->dft->filter->idft-vindue
-Bjarne

---

On Tue, 17 Apr 2001 00:10:48 +0200, Bjarne Laursen wrote:
>Nej det er nok ikke nogen god ide. Hvis du kun har en blok skal du
>sørge for at signalet kan være i blokken uden at nå endepunkterne, så
>behøver du ingen extra vinduefunktion og så kan signalet gendannes med
>idft.
>Hvis du derimod har delt dit signal op i flere blokke (sådan som
>diverse audiofiltre virker) så skal du så vidt jeg ved lade blokkene
>overlappe. Vinduefunktionen skal så være således at når du summere
>blokkene bagefter giver det 1. Det er sikkert en god ide at bruge
>vinduefunktionen både før og efter: vindue->dft->filter->idft-vindue

Mit signal er ikke opdelt. Det er en ganske almindelig DFT
og IDFT (skulle jeg mene), men det virker ikke. Er det
rigtigt, at jeg skal arbejde med komplekse talpar? Skal det
i stedet bare være reelle?

Jeg må godt nok sige, at der ikke står særligt meget om
implementering af dette i mine bøger.. =\.

http://stephan.tisprut.dk/Skrammel/fourier/fourier.cpp



--
Stephan Henningsen /
/ http://tisprut.dk

---

stephan@sanxion.tisprut.dk (Stephan Henningsen) wrote:

>Er det
>rigtigt, at jeg skal arbejde med komplekse talpar? Skal det
>i stedet bare være reelle?
Et array med complexe tal skulle være ok.
Jeg har kigget på din kode og jeg synes ikke det ser helt rigtigt ud:

>         c.re += (X[k].re + X[k].im) * cos( (2*PI/N)*k*n );
>         c.im += (X[k].re + X[k].im) * sin( (2*PI/N)*k*n );
Så vidt jeg husker skal du multiplicere to complex tal. Det ser
således ud:

c=a*b
<=>
c.re = a.re*b.re - a.im*b.im;
c.im= a.re*b.im + a.im*b.re;

Altså:

c.re += X[k].re*cos( (2*PI/N)*k*n ) - X[k].im*sin( (2*PI/N)*k*n );
c.im += X[k].re*sin( (2*PI/N)*k*n ) + X[k].im*cos( (2*PI/N)*k*n );


Hvis dine X[k].im altid er nul så får du det samme som før, men det er
de nok næppe i din idft.

-Bjarne

---

On Tue, 17 Apr 2001 01:37:48 +0200, Bjarne Laursen wrote:
>Altså:
>
>c.re += X[k].re*cos( (2*PI/N)*k*n ) - X[k].im*sin( (2*PI/N)*k*n );
>c.im += X[k].re*sin( (2*PI/N)*k*n ) + X[k].im*cos( (2*PI/N)*k*n );

Du har ret. Jeg fandt også ud af, at N ikke skal køre til
N/2, men til hele N. De N/2 er kun, hvis man ikke ønsker at
genskabe signalet, da den sidste halvdel er mere eller
mindre redundant. Men for at en IDFT skal virke, skal jeg
have hele signalet.

Nu virker min IDFT da også -- men kun på første bufferplads, X[0];
resten af bufferen bliver transformeret til noget mærkeligt
med meget, meget små værdier:

Read 10 samples from file test1.dec successfully.
x[0] = ( 0; 0)
x[1] = ( 1; 0)
x[2] = ( 2; 0)
x[3] = ( 3; 0)
x[4] = ( 4; 0)
x[5] = ( 5; 0)
x[6] = ( 6; 0)
x[7] = ( 7; 0)
x[8] = ( 8; 0)
x[9] = ( 9; 0)
Fourier-transforming: 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% done.
X_dft[0] = ( 45; 0)
X_dft[1] = ( -5;15.388418)
X_dft[2] = ( -5;6.8819096)
X_dft[3] = ( -5;3.6327126)
X_dft[4] = ( -5;1.6245985)
X_dft[5] = ( -5;-5.5107286e-15)
X_dft[6] = ( -5;-1.6245985)
X_dft[7] = ( -5;-3.6327126)
X_dft[8] = ( -5;-6.8819096)
X_dft[9] = ( -5;-15.388418)
Inverse Fourier-transforming: 0% 10% 20% 30% 40% 50% 60% 70% 80% 90% done.
x_idft[0] = ( 0; 0)
x_idft[1] = (-1.4815623e-16;6.1203213e-17)
x_idft[2] = (-8.5516446e-16;5.1770654e-16)
x_idft[3] = (-2.4318655e-16;1.1158609e-15)
x_idft[4] = (-8.6779542e-16;2.3019781e-15)
x_idft[5] = ( 0;3.0615159e-15)
x_idft[6] = (1.5087757e-15;4.6126297e-15)
x_idft[7] = (4.1212693e-15;6.0069137e-15)
x_idft[8] = (1.057379e-14;6.9380265e-15)
x_idft[9] = (3.0428784e-14;9.924353e-15)
Wrote 10 samples to file buffer.idft successfully.

Mystisk. Min IDFT og DFT kan ses på
http://stephan.tisprut.dk/Skrammel/fourier/fourier.cpp





--
Stephan Henningsen /
/ http://tisprut.dk

---

On Mon, 16 Apr 2001 16:08:11 +0000, stephan@sanxion.tisprut.dk
(Stephan Henningsen) wrote:


>>>Men jeg problemer med min IDFT. Skal en invers
>>>F-transformation ikke genskabe signalet, som det så ud, før
>>>det blev F-transformeret? Skal man vægte det med en invers
>>>Hanning?
>>
>int
>fourier :: idft ( const vector<complex>& X, vector<complex>& x )
>{
> unsigned long int k, K, n, N;
> complex c;
>
>
> x.clear(); // Clear any existing vector before pushing data into it.
> K = X.size();
> N = K*2;
> for ( n = 0 ; n < N ; n++ ) {
> c.re = c.im = 0;
>
> for ( k = 0 ; k < N ; k++ ) {
> c.re += (X[k].re + X[k].im) * cos( (2*PI/N)*k*n );
> c.im += (X[k].re + X[k].im) * sin( (2*PI/N)*k*n );
> }
>
> c.re /= N;
> c.im /= N;
>
> x.push_back( c );
> }
>
>
> return 0;
>}

Hej igen

Jeg har kigget lidt på din IDFT, og fundet lidt fejl i din
udregning.....

c.re += X[k].re*cos(2*PI/N*k*n) - X[k].im*sin(2*PI/N*k*n);
c.im += X[k].re*sin(2*PI/N*k*n) + X[k].im*cos(2*PI/N*k*n);

Dette er rigtigt... Iøvrigt skal du for hvert et k, udregne alle N
punkter.

Det har aldrig været et problem for mig.

Her en lille opfriskning af DFT og IDFT...

DFT:
X[k] = sum( 0, N-1, x[n]*exp(-j*2*PI/N*k*n) )


IDFT:
x[n] = 1/N * sum( 0, N-1, X[k]*exp(j*2*PI/N*k*n) )

Husk at x[n] og X[k] kan begge være komplekse..

/Thomas

---

On Sat, 21 Apr 2001 12:46:46 +0200, Thomas Lykkeberg wrote:
>Jeg har kigget lidt på din IDFT, og fundet lidt fejl i din
>udregning.....
>
>c.re += X[k].re*cos(2*PI/N*k*n) - X[k].im*sin(2*PI/N*k*n);
>c.im += X[k].re*sin(2*PI/N*k*n) + X[k].im*cos(2*PI/N*k*n);

Jeps, jeg har selv igen prøvet at regne på det efter en
opfriskning i de komplekse regneregler, og kom frem til
netop det samme.


>Dette er rigtigt... Iøvrigt skal du for hvert et k, udregne alle N
>punkter.

Ja, det fandt jeg ud af efter en snak med min underviser.


>Det har aldrig været et problem for mig.

Jeg fatter heller ikke, hvad der går galt her, og jeg er ved
at være seriøst træt af det. .. Hov, nu virker det vist.


>Her en lille opfriskning af DFT og IDFT...
>
>DFT:
>X[k] = sum( 0, N-1, x[n]*exp(-j*2*PI/N*k*n) )
>
>
>IDFT:
>x[n] = 1/N * sum( 0, N-1, X[k]*exp(j*2*PI/N*k*n) )
>
>Husk at x[n] og X[k] kan begge være komplekse..

Jeps, den er jeg fuldstændig med på. Jeg vil lige piske
nogle signaler igennem mit program og se om det virker, som
det skal.. jeg vender frygteligt tilbage =)


--
Stephan Henningsen /
/ http://tisprut.dk

---

On Sat, 21 Apr 2001 12:46:46 +0200, Thomas Lykkeberg wrote:
>Hej igen
>
>Jeg har kigget lidt på din IDFT, og fundet lidt fejl i din
>udregning.....
>
>c.re += X[k].re*cos(2*PI/N*k*n) - X[k].im*sin(2*PI/N*k*n);
>c.im += X[k].re*sin(2*PI/N*k*n) + X[k].im*cos(2*PI/N*k*n);

Jeg kørte det for X[n] og ikke X[k]. Sådan kan det gå, når
man er for hurtig med cut-n-paste =(.

Nu virker det perfekt. Tak for hjælpen!

Så skal der lige skrives et program, der kan lave nogle
wav-filer udfra en buffer =).



--
Stephan Henningsen /
/ http://tisprut.dk

---

Stephan Henningsen wrote:

> Så skal der lige skrives et program, der kan lave nogle
> wav-filer udfra en buffer =).

Led efter "sox" på nettet, og overvej at give den dine rå data filer.

--
Thorbjørn Ravn Andersen "...plus... Tubular Bells!"
http://bigfoot.com/~thunderbear

---

On Tue, 24 Apr 2001 07:31:00 +0200, Thorbjørn Ravn Andersen wrote:
>Stephan Henningsen wrote:
>
>> Så skal der lige skrives et program, der kan lave nogle
>> wav-filer udfra en buffer =).
>
>Led efter "sox" på nettet, og overvej at give den dine rå data filer.

Jeps, der er styr på det. Jeg havde i starten nogle
problemer med signed og unsigned, men det spiller nu. Men
min write_raw_file-funktioner får filerne til at fylde to
bytes mere (unsigned short int), altså én sample mere. Data
i starten og slutningen den oprindelige fil og kopien er
ens, så der må gå galt et sted midt i filen. Jeg må lige se
lidt nærmere på det. Men det er mystisk.


--
Stephan Henningsen /
/ http://tisprut.dk

---

Den Sun, 29 Apr 2001 02:02:15 +0000 skrev Stephan Henningsen:
>On Tue, 24 Apr 2001 07:31:00 +0200, Thorbjørn Ravn Andersen wrote:
>>Stephan Henningsen wrote:
>>
>>> Så skal der lige skrives et program, der kan lave nogle
>>> wav-filer udfra en buffer =).
>>
>>Led efter "sox" på nettet, og overvej at give den dine rå data filer.
>
>Jeps, der er styr på det. Jeg havde i starten nogle
>problemer med signed og unsigned, men det spiller nu. Men
>min write_raw_file-funktioner får filerne til at fylde to
>bytes mere (unsigned short int), altså én sample mere. Data
>i starten og slutningen den oprindelige fil og kopien er
>ens, så der må gå galt et sted midt i filen. Jeg må lige se
>lidt nærmere på det. Men det er mystisk.

Prøv at køre "cmp fil1 fil2" - den burde fortælle hvor forskellen er.

Mvh
Kent
--
Nu med en e-mail adresse der virker...

---

---

On Mon, 30 Apr 2001 15:14:01 +0200, Thorbjoern Ravn Andersen wrote:
>On Sun, 29 Apr 2001, Stephan Henningsen wrote:
>
>Har du husket et "b" når du åbner for skrivning?
>
>I en sådan situation vil jeg under Unix køre "od -c" på begge filer og så
>køre diff på det. Første linie i uddata er hvor de begynder at være
>forskellige..

Hmm, der er ganske rigtigt en forskel i slutningen af den
fil, jeg selv skriver.

Hvad mener du med det "b"? Er dette ikke godt nok:


ofstream f(filename, ios::binary);

   [...]

N = x.size();

   [...]

if ( !f.is_open() ) {
cerr << "write_raw_buffer: unable to write raw file: "
<< filename << endl;
return -1;
}

for ( n = 0 ; n < N ; n++ ) {
i = (unsigned short int) rint( x[n].re );
f.write( &i, sizeof(i) );
}
f.close();

--
Stephan Henningsen /
/ http://tisprut.dk

---

Stephan Henningsen wrote:

> Hmm, der er ganske rigtigt en forskel i slutningen af den
> fil, jeg selv skriver.
>
> Hvad mener du med det "b"? Er dette ikke godt nok:

Aner jeg ikke - jeg er ikke C++ programmør.

I rå C benytter man "b" som parameter under non-Unix miljøer til at
fortælle at man IKKE vil have gjort noget som helst ved de bytes der
skal skrives.

Hvad præcis er forskellen.
--
Thorbjørn Ravn Andersen "...plus... Tubular Bells!"
http://bigfoot.com/~thunderbear

---

"Stephan Henningsen" <stephan@sanxion.tisprut.dk> wrote in message
news:slrn9d9lpi.9i.stephan@sanxion.tisprut.dk...
> On Wed, 11 Apr 2001 18:52:58 +0200, Thomas Lykkeberg wrote:
> >
> >Der sker nemlig det når du sampler et signal, at du blot omsætter den
> >samplede "stump" til et laaaaangt repeterende signal. Jeg ved ikke om
> >du kan følge mig her.?
>
> Jeg er nogenlunde med. Så for at kompenerer for, at det
> ikke er et repeterende signal man sampler, men et endeligt,
> hvilket giver visser fejl i udregningen, så folder man en
> slags nødløsning på?

Ja. DFT kræver at signalet er gentager sig, både i frekvens og tid. Da
en fysisk måling ikke klan forventes at være et præcist udsnit, så
ganger man et vindue på. Gø man ikke det, så kan der være diskontuinitet
ved kanterne. Dette vil følge til overtoner, lidt ligesom
firkant-signaler består af en masse overtoner.


--
Morten Boysen