图象傅立叶转换.txt

一个图象傅立叶转换

#include <conio.h>
#include <stdio.h>
#include <stdlib.h>
#include <process.h>
#include <dos.h>
#include <dir.h>
#include <io.h>
#include <string.h>
#include <fcntl.h>
#include <alloc.h>
#include <math.h>
#include <graphics.h>

union REGS       iregs, oregs, sregs;
double data1[512], data2[512], max1, bu1[3], bu2[3];
unsigned char image1[128][128], image2[128][128];
int max;

/* isi=-1   Fourier transfer;
   isi=1    Reverse Fourier transfer;
   n:      Total point number is 2* *m.
*/
void  fft ( int n, int isi )
{
        int nn, m, i, j, mmax, istep;
        double pi, fn, temp1, temp2, theta, w1, w2;
        nn=1;

        for(i=1;i<15;i++){
                m=i;
                nn=nn*2;
                if(nn==n)break;
        }
        pi=4.0*atan(1.0);
        fn=(double)n;
        j=1;
/* Reverse bit */
        for(i=1;i<=n;i++){
                if((i-j)<0){
                        temp1=data1[j-1];
                        temp2=data2[j-1];
                        data1[j-1]=data1[i-1];
                        data2[j-1]=data2[i-1];
                        data1[i-1]=temp1;
                        data2[i-1]=temp2;
                }
                m=n/2;
ld50:           if((j-m)>0){
                        j=j-m;
                        m=(m+1)/2;
                        goto ld50;
                }
                j=j+m;
        }
        mmax=1;
ld90:   if(mmax-n<0){
                istep=2*mmax;
                for(m=1;m<=mmax;m++){
                        theta=pi*(double)(isi*(m-1))/(double)mmax;
                        w1=cos(theta);
                        w2=sin(theta);
                        for(i=m;i<n;i+=istep){
                                j=i+mmax;
                                temp1=w1*data1[j-1]-w2*data2[j-1];
                                temp2=w1*data2[j-1]+w2*data1[j-1];
                                data1[j-1]=data1[i-1]-temp1;
                                data2[j-1]=data2[i-1]-temp2;
                                data1[i-1]=data1[i-1]+temp1;
                                data2[i-1]=data2[i-1]+temp2;
                        }
                }
                mmax=istep;
                goto ld90;
        }
        if(isi>=0){
                for(i=0;i<n;i++){
                        data1[i]=data1[i]/fn;
                        data2[i]=data2[i]/fn;
                }
        }
}

main()
{
	FILE *fp, *fp1;

        int a,b,c,x,y,z;
        int n, isi, driver, mode1;
	unsigned char d[128][128];
        driver=DETECT;
        mode1=2;
        initgraph(&driver,&mode1,"");

        iregs.x.ax=0x0100;
        iregs.x.cx=0x0000;
        int86(0x10,&iregs,&oregs);

        closegraph();

        driver=DETECT;
        mode1=2;
        initgraph(&driver,&mode1,"");

        iregs.x.ax=0x005d;
        int86(0x10,&iregs,&oregs);

        for(a=0;a<256;a++){
            iregs.x.ax=0x1010;
            iregs.x.bx=a;
            iregs.h.dh=a/4;
            iregs.h.ch=a/4;
            iregs.h.cl=a/4;
            int86(0x10,&iregs,&oregs);
        }

	fp=fopen("buffer","wb+");
	fp1=fopen("buffer1","wb+");
        for(a=0;a<128;a++){
                for(b=0;b<128;b++) {
			if(a>=40&&a<88&&b>=40&&b<88) data1[b]=100.0;
			else data1[b]=0.0;
			data2[b]=0.0;
			putpixel2(b,a+200,255-(int)data1[b]);
                }
                fft( 128, -1 );
		fwrite(data1,sizeof(double),128,fp);
                fwrite(data2,sizeof(double),128,fp);
	}
	fclose(fp);

	fp=fopen("buffer","rb");
	max=0;
	for(a=0;a<128;a++){
		for(b=0;b<128;b++){
			fseek(fp,(long)b*2048+a*8,SEEK_SET);
			fread(bu1,sizeof(double),1,fp); data1[b]=bu1[0];
			fseek(fp,(long)b*2048+a*8+1024,SEEK_SET);
			fread(bu2,sizeof(double),1,fp); data2[b]=bu2[0]; 
                }
		fft( 128, -1 );

		for(b=0;b<128;b++){
			data1[b]=sqrt(data1[b]*data1[b]+data2[b]*data2[b]);
			data1[b]=log(data1[b]+1.0);
		}

		for(b=0;b<128;b++) if(max<data1[b])max=data1[b];
		fwrite(data1,sizeof(double),128,fp1);
	}
	fclose(fp);
	fclose(fp1);

	fp=fopen("buffer1","rb");
	fp1=fopen("buf","wb");
	for(a=0;a<128;a++){
		fread(data1,sizeof(double),128,fp);
		for(b=0;b<128;b++){
		     data1[b]=data1[b]*256.0/(max+1);
		     putpixel2((a+64)%128,(b+64)%128,(int)data1[b]);
		     d[(a+64)%128][(b+64)%128]=(unsigned char)data1[b]; 
		 }
	 }
	 for(a=0;a<128;a++)
		fwrite(d[a],sizeof(unsigned char),128,fp1);
	fclose(fp);
	fclose(fp1);
        getch();
	closegraph();

}

getpixel2( int a, int b )
{
    iregs.h.ah=0xd;
    iregs.x.dx=b;
    iregs.x.cx=a;
    iregs.x.bx=0;
    int86(0x10,&iregs,&oregs);
    return(oregs.h.al);
}

putpixel2( int a, int b,int num)
{
    iregs.h.ah=0xc;
    iregs.h.al=num;
    iregs.x.dx=b;
    iregs.x.cx=a;
    iregs.x.bx=0;
    int86(0x10,&iregs,&oregs);
}