<![CDATA[#include <stdio.h>
#include <errno.h>
#include <math.h>
#define N 1024L
logtbl[] = { 0, 30, 48, 60, 70, 78, 85, 90, 96, 100, 104, 108,
112, 115, 118, 121, 123, 126, 128, 130, 132, 134, 136, 138,
140, 142, 143, 145, 146, 148, 149, 151, 152, 153, 155, 156,
157, 158, 159, 160, 161, 163, 164, 165, 165, 166, 167, 168,
169, 170, 171, 172, 173, 173, 174, 175, 176, 177, 177, 178,
179, 179, 180, 181, 181, 182, 183, 183, 184, 185, 185, 186,
187, 187, 188, 188, 189, 189, 190, 190, 191, 192, 192, 193,
193, 194, 194, 195, 195, 196, 196, 197, 197, 197, 198, 198,
199, 199, 200, 200, 201, 201, 201, 202, 202, 203, 203, 204,
204, 204, 205, 205, 205, 206, 206, 207, 207, 207, 208, 208,
208, 209, 209, 210, 210, 210, 211, 211, 211, 212, 212, 212,
213, 213, 213, 214, 214, 214, 214, 215, 215, 215, 216, 216,
216, 217, 217, 217, 217, 218, 218, 218, 219, 219, 219, 219,
220, 220, 220, 221, 221, 221, 221, 222, 222, 222, 222, 223,
223, 223, 223, 224, 224, 224, 224, 225, 225, 225, 225, 226,
226, 226, 226, 227, 227, 227, 227, 228, 228, 228, 228, 229,
229, 229, 229, 229, 230, 230, 230, 230, 230, 231, 231, 231,
231, 232, 232, 232, 232, 232, 233, 233, 233, 233, 233, 234,
234, 234, 234, 234, 235, 235, 235, 235, 235, 236, 236, 236,
236, 236, 237, 237, 237, 237, 237, 237, 238, 238, 238, 238,
238, 239, 239, 239, 239, 239, 239, 240, 240, 240, 240, 240,
240, 241, 241
};
main()
{
FILE *fp;
int isclf, i;
short *idata;
short *out;
idata = (short *)malloc(2*N*sizeof(short));
if (idata==-1) {
fprintf(stderr, "No memory\n");
exit(errno);
}
out = (short *)malloc(2*N*sizeof(short));
if (out==-1) {
fprintf(stderr, "No memory\n");
exit(errno);
}
fp = fopen("/n0/eddy/dd/fft/spectrum.dat", "r");
if (!fp) {
fprintf(stderr, "Can't open spectrum.dat\n");
exit(errno);
}
fread(idata, sizeof(short), 2*N, fp);
initdisplay();
isclf = ffti2(N, 0L, idata, &idata[N]);
if (isclf!=-1) {
normalize(idata, &idata[N], out, isclf);
display(out, isclf);
}
else
printf("error in ffti2\n");
}
initdisplay()
{
int i;
cls(0);
for (i=0; i<512; i+=51)
line(i,10,i,255,20);
for (i=10; i<256; i+=24)
line(0,i,511,i,20);
}
display(out, isclf)
short *out;
int isclf;
{
short i, lp, y1, y2;
int k = 6+isclf;
y1 = 0;
for (i=0; i<N/2; i++) {
lp = (lp+*out++>>k)>>1;
y2 = logtbl[lp];
vline(i, y1, y2);
y1 = y2;
}
}
vline(i, y1, y2)
int i, y1, y2;
{
int j;
if (y2>=y1)
for(j=y1; j<=y2; j++)
xpoint(i, j);
else
for (j=y2; j<y1; j++)
xpoint(i, j);
}
normalize(re, im, out)
short *re, *im, *out;
{
int i;
register short x, y;
*out++ = re[0];
for (i=1; i<N/2; i++) {
x = (re[i]+re[N-i])>>1;
y = (im[i]-im[N-i])>>1;
*out++ = fsqrt(x*x+y*y);
}
*out++ = im[0];
for (i=1; i<N/2; i++) {
x = (im[i]+im[N-i])>>1;
y = (re[N-i]-re[i])>>1;
*out++ = fsqrt(x*x+y*y);
}
}
#asm
* bereken wortel
fsqrt: movem.l d2-d4,-(sp)
move.l d0,d2
beq w30
cmp.l #1,d2
beq w30
cmp.l #$ffff,d2
ble w10
move #$7fff,d0
bra w20
w10 move d2,d0
lsr #1,d0
w20 move d0,d3
move.l d2,d4
divu d0,d4
add d4,d0
lsr #1,d0
cmp d0,d3
bgt w20
w30 movem.l (sp)+,d2-d4
rts
#endasm
(68020) Super: ]]>