unit FFT;

interface

	uses
		QuickDraw, Palettes, PrintTraps, globals, Utilities;
	const
		NNNx2 = 1024;
		NNNN = 150000; { 256 * 290 * 2 }
{ so krumm wegen rainer }



	type
		FFTVector = array[0..NNNx2] of real;
		FFTArray = array[0..NNNN] of real;
		FFTArrayPtr = ^FFTArray;
		FFTArrayHandle = ^FFTArrayPtr;
		CorrBut = (AutoB, Fixed);
		FFT_Info = record
				FFT_Done: boolean;
				FFTPict_Info: InfoPtr;
			end;
		SavedPict = array[0..NNNN] of integer;
		SavedPictPtr = ^SavedPict;
		SavedPictHandle = ^SavedPictPtr;



	function DuplicateFFT (SavingBlankField: boolean; var w, h: longint): boolean;
	function DuplicateHorizontalFFT (SavingBlankField: boolean; var w, h: longint): boolean;
	procedure DoHorizontalFFT;
	procedure DoHorizontalBackFFT;
	procedure Do2DFFT;
	procedure Do2DBackFFT;
	procedure UndoFFT;
	procedure Restore;
	procedure DoFFTFilter (n: integer);
	procedure DoBusyBox (n, a: integer);
	procedure OpenBusyBox (text: str255);
	procedure CloseBusyBox;
	procedure SmoothPicture (FFTArrayH, ReFFTArrayH: FFTArrayHandle; n: longint);
	procedure ShowFFTPhases;
	procedure ReadPictandScale (FFTArrayH: FFTArrayHandle; n: longint);
	procedure ReadPictHorizontal (FFTArrayH: FFTArrayHandle; n, height: longint);
	procedure ReadPictAndScaleHorizontal (FFTArrayH: FFTArrayHandle; n, height: longint);
	procedure FFTHorizontal (FFTArrayh: FFTArrayhandle; n, height, isign: longint);
	procedure WriteHorizontalFFTPict (FFTArrayH: FFTArrayHandle; n, height: longint);
	procedure CopyFFTArray (FFTArrayH, ReFFTArrayH: FFTArrayHandle; w, h: longint);
	procedure ReFFTHorizontal (ReFFTArrayH: FFTArrayhandle; w, h, isign: longint);
	procedure WriteFFTPict (FFTArrayH: FFTArrayHandle; n: longint);
	procedure WriteRealPict (ReFFTArrayH: FFTArrayHandle; w, h: longint);
	procedure WriteAndScaleRealPict2D (ReFFTArrayH: FFTArrayHandle; w, h: longint; ScalingFactor: real; Button: CorrBut);
	procedure WriteAndScaleRealPictHorizontal (ReFFTArrayH: FFTArrayHandle; w, h: longint; ScalingFactor: real; Button: CorrBut);
	procedure FFT2D (FFTArrayH: FFTArrayHandle; n, isign: longint);
	procedure ReFFT2D (ReFFTArrayH: FFTArrayHandle; n, isign: longint);

	const
		NNN = 256;
		BusyDLOG = 2200;
		BusyBox = 3;

	var
		FFTArrayH: FFTArrayHandle;
		ReFFTArrayH: FFTArrayHandle;
		FFTInfo: FFT_Info;
		SavedPictH: SavedPictHandle;
		FFTFilterNoise, FFTFilterComponents: integer;


implementation

	const
		FilterDLOG = 2100;
		FilterOk = 1;
		FilterCancel = 2;
		NoisePercent = 5;
		NoiseUp = 6;
		NoiseDown = 7;
		FlatenText = 10;
		FlatenUp = 11;
		FlatenDown = 12;
		NoiseFastUp = 13;
		NoiseFastDown = 14;
		FlatenFastUp = 15;
		FlatenFastDown = 16;

	var
		busyDialog: DialogPtr;
		BusyRect: Rect;


	function DuplicateFFT;{(SavingBlankField:boolean)}
		var
			name: str255;
			hstart, vstart: integer;
			i: integer;
			SaveInfo: InfoPtr;
			src, dst: ptr;
			offset: longint;
			AutoSelectAll: boolean;
			Size: longint;
			isOkay: boolean;
	begin
		isOkay := false;
		WhatToUndo := NothingToUndo;
		if (not SavingBlankField) and (NotRectangular or NotinBounds) then begin
				DuplicateFFT := false;
				exit(DuplicateFFT);
			end;
		AutoSelectAll := (not Info^.RoiShowing) or SavingBlankField;
		if AutoSelectAll then
			SelectAll(false);
		ShowWatch;
		with info^ do begin
				if SavingBlankField then
					name := 'Blank Field'
				else begin
						GetWTitle(wptr, name);
						if pos('FFT - ', name) = 1 then
							name := concat('Re', name);
						if pos('FFT - ', name) = 0 then
							name := concat('FFT - ', name);
						if length(name) > 32 then
							delete(name, 33, length(name) - 32);
					end;
				with RoiRect do begin
						w := right - left;
						if odd(w) and (left + w < PicRect.right) then
							w := w + 1;
						h := bottom - top;
						hstart := left;
						vstart := top;
					end;
				if h > w then
					h := w
				else
					w := h;
				begin
					Size := round(h div 32);
					case Size of
						0, 1: 
							h := 32;
						2, 3: 
							h := 32 * 2;
						4, 5, 6, 7: 
							h := 32 * 4;
						8, 9, 10, 11, 12, 13, 14, 15: 
							h := 32 * 8;
{16: }
{h := 32 * 16;}
						otherwise
							h := NNN;
					end;
				end;
				w := h;
			end;
		if AutoSelectAll then
			KillRoi;
		SaveInfo := Info;
		if NewPicWindow(name, w, h) then
			with SaveInfo^ do begin
					Info^.x_range := x_range;
					Info^.y_range := y_range;
					Info^.z_range := z_range;

					offset := longint(vstart) * BytesPerRow + hstart;
					src := ptr(ord4(PicBaseAddr) + offset);
					dst := Info^.PicBaseAddr;
					for i := 0 to h - 1 do begin
							BlockMove(src, dst, w);
							src := ptr(ord4(src) + BytesPerRow);
							dst := ptr(ord4(dst) + w);
						end;
					if SavingBlankField then begin
							Info^.PIctureType := BlankField;
							BlankFieldInfo := info;
						end;
					isOkay := true;
				end
		else
			isOkay := false;
		with Info^ do
			ScaleToFitWindow := true;
		DuplicateFFT := isOkay;
	end;

	function DuplicateHorizontalFFT;{(SavingBlankField:boolean)}
		var
			name: str255;
			width, height, hstart, vstart: integer;
			i: integer;
			SaveInfo: InfoPtr;
			src, dst: ptr;
			offset: longint;
			AutoSelectAll: boolean;
			Size: longint;
			isOkay: boolean;
	begin
		isOkay := false;
		WhatToUndo := NothingToUndo;
		if (not SavingBlankField) and (NotRectangular or NotinBounds) then begin
				duplicateHorizontalFFT := false;
				exit(DuplicateHorizontalFFT);
			end;
		AutoSelectAll := (not Info^.RoiShowing) or SavingBlankField;
		if AutoSelectAll then
			SelectAll(false);
		ShowWatch;
		with info^ do begin
				if SavingBlankField then
					name := 'Blank Field'
				else begin
						GetWTitle(wptr, name);
						if pos('FFT - ', name) = 1 then
							name := concat('Re', name);
						if pos('FFT - ', name) = 0 then
							name := concat('FFT - ', name);
						if length(name) > 32 then
							delete(name, 33, length(name) - 32);
					end;
				with RoiRect do begin
						width := right - left;
						if odd(width) and (left + width < PicRect.right) then
							width := Width + 1;
						height := bottom - top;
						if odd(height) then
							height := height - 1;
						hstart := left;
						vstart := top;
					end;
				begin
					Size := round(width div 32);
					case Size of
						0, 1: 
							width := 32;
						2, 3: 
							width := 32 * 2;
						4, 5, 6, 7: 
							width := 32 * 4;
						8, 9, 10, 11, 12, 13, 14, 15: 
							width := 32 * 8;
						16, 17, 18, 19, 20: 
							width := 32 * 16;
						otherwise
							width := NNN;
					end;
				end;
				w := width;
				h := height;
			end;
		if AutoSelectAll then
			KillRoi;
		SaveInfo := Info;
		if NewPicWindow(name, width, height) then
			with SaveInfo^ do begin
					Info^.x_range := x_range;
					Info^.y_range := y_range;
					Info^.z_range := z_range;

					offset := longint(vstart) * BytesPerRow + hstart;
					src := ptr(ord4(PicBaseAddr) + offset);
					dst := Info^.PicBaseAddr;
					for i := 0 to height - 1 do begin
							BlockMove(src, dst, width);
							src := ptr(ord4(src) + BytesPerRow);
							dst := ptr(ord4(dst) + width);
						end;
					if SavingBlankField then begin
							Info^.PIctureType := BlankField;
							BlankFieldInfo := info;
						end;
					isOkay := true;
				end
		else
			isOkay := false;
		with Info^ do
			ScaleToFitWindow := true;
		DuplicateHorizontalFFT := isOkay
	end;


	procedure BitRev (var Dat: FFTVector; n: longint);
		var
			i, j, k, m, n2: longint;
			tr, ti: real;
	begin
		n2 := 2 * n;
		j := 0;
		for i := 0 to n2 - 1 do begin
				k := 2 * i;
				if (j > k) then begin
						tr := Dat[j];
						ti := Dat[j + 1];
						Dat[j] := Dat[k];
						Dat[j + 1] := Dat[k + 1];
						Dat[k] := tr;
						Dat[k + 1] := ti;
					end;
				m := n2 div 2;
				while ((m >= 2) and (j >= m)) do begin
						j := j - m;
						m := m div 2;
					end;
				j := j + m;
			end;
	end;



	procedure FFT (var Dat: FFTVector; n, isign: longint);
		var
			i, is, j, m, mmax, n2: longint;
			tr, ti, wr, wpr, wi, wpi, wt, theta, aux: real;
	begin
		BitRev(Dat, n);
		n2 := 2 * n;
		mmax := 2;
		while (n2 > mmax) do begin
				is := 2 * mmax;
				theta := (2 * Pi) / (isign * mmax);
				wpr := -2.0 * sqr(sin(0.5 * theta));
				wpi := sin(theta);
				wr := 1.0;
				wi := 0.0;
				m := 0;
				while (m < mmax) do begin
						i := m;
						while (i < n2) do begin
								j := i + mmax;
								tr := wr * Dat[j] - wi * Dat[j + 1];
								ti := wr * Dat[j + 1] + wi * Dat[j];
								Dat[j] := Dat[i] - tr;
								Dat[j + 1] := Dat[i + 1] - ti;
								Dat[i] := Dat[i] + tr;
								Dat[i + 1] := Dat[i + 1] + ti;
								i := i + is;
							end;
						wt := wr;
						wr := wr * wpr - wi * wpi + wr;
						wi := wi * wpr + wt * wpi + wi;
						m := m + 2;
					end;
				mmax := is;
			end;
		if (isign = 1) then
			for i := 0 to n2 - 1 do
				Dat[i] := Dat[i] / n;
	end;



	procedure twofft (var fft1, fft2: FFTVector; n, isign: longint);
 (* Programs using routine TWOFFT must define types}
{TYPE}
{   FFTVector = ARRAY [0..2*n-1] OF real;}
{where n is the dimension of the real-valued data arrays. *)
		var
			nn3, nn2, nn, jj, j: integer;
			rep, rem, aip, aim: real;
	begin
		nn := n + n;
		nn2 := nn + 2;
		nn3 := nn + 3;
		FFT(fft1, n, isign);
		fft2[0] := fft1[1];
		fft1[1] := 0.0;
		fft2[1] := 0.0;
		for jj := 1 to (n div 2) do begin
				j := 2 * jj + 1;
				rep := 0.5 * (fft1[j - 1] + fft1[nn2 - j - 1]);
				rem := 0.5 * (fft1[j - 1] - fft1[nn2 - j - 1]);
				aip := 0.5 * (fft1[j] + fft1[nn3 - j - 1]);
				aim := 0.5 * (fft1[j] - fft1[nn3 - j - 1]);
				fft1[j - 1] := rep;
				fft1[j] := aim;
				fft1[nn2 - j - 1] := rep;
				fft1[nn3 - j - 1] := -aim;
				fft2[j - 1] := aip;
				fft2[j] := -rem;
				fft2[nn2 - j - 1] := aip;
				fft2[nn3 - j - 1] := rem
			end
	end;

	procedure FFTHorizontal (FFTArrayH: FFTArrayhandle; n, height, isign: longint);
		var
			i, j: longint;
			LineIndex, CollumnIndex: longint;
			nx2, nx4: longint;
			n_1: longint;
			n_half: longint;
			FFTDataPtr: FFTArrayPtr;
			fft1, fft2: FFTVector;
	begin
		OpenBusyBox('Processing Horizontal Fourier Transform');
		nx2 := n * 2;
		nx4 := n * 4;
		n_1 := n - 1;
		n_half := trunc(height / 2);
		LineIndex := 0;
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		for i := 0 to n_half - 1 do begin
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						fft1[CollumnIndex] := FFTDataPtr^[LineIndex + CollumnIndex];
						fft1[CollumnIndex + 1] := FFTDataPtr^[LineIndex + CollumnIndex + nx2];
						CollumnIndex := CollumnIndex + 2;
					end;
				twofft(fft1, fft2, n, isign);
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						FFTDataPtr^[LineIndex + CollumnIndex] := fft1[CollumnIndex];
						FFTDataPtr^[LineIndex + CollumnIndex + 1] := fft1[CollumnIndex + 1];
						FFTDataPtr^[LineIndex + CollumnIndex + nx2] := fft2[CollumnIndex];
						FFTDataPtr^[LineIndex + CollumnIndex + nx2 + 1] := fft2[CollumnIndex + 1];
						CollumnIndex := CollumnIndex + 2;
					end;
				DoBusyBox(i * 2, height);
				LineIndex := LineIndex + nx4;
			end;
		HUnlock(handle(FFTArrayH));
		CloseBusyBox;
	end;

	procedure FFT2D (FFTArrayH: FFTArrayhandle; n, isign: longint);
		var
			i, j: longint;
			LineIndex, CollumnIndex: longint;
			nx2: longint;
			nx4: longint;
			DatVector: FFTVector;
			FFTDataPtr: FFTArrayPtr;
			nxnx2: longint;
			n_1: longint;
			n_half: longint;
			fft1, fft2: FFTVector;
	begin
		OpenBusyBox('Processing Fourier Transform');
		;
		nx2 := n * 2;
		nx4 := n * 4;
		nxnx2 := n * n * 2;
		n_1 := n - 1;
		n_half := trunc(n / 2);
		LineIndex := 0;
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		for i := 0 to n_half - 1 do begin
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						fft1[CollumnIndex] := FFTDataPtr^[LineIndex + CollumnIndex];
						fft1[CollumnIndex + 1] := FFTDataPtr^[LineIndex + CollumnIndex + nx2];
						CollumnIndex := CollumnIndex + 2;
					end;
				twofft(fft1, fft2, n, isign);
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						FFTDataPtr^[LineIndex + CollumnIndex] := fft1[CollumnIndex];
						FFTDataPtr^[LineIndex + CollumnIndex + 1] := fft1[CollumnIndex + 1];
						FFTDataPtr^[LineIndex + CollumnIndex + nx2] := fft2[CollumnIndex];
						FFTDataPtr^[LineIndex + CollumnIndex + nx2 + 1] := fft2[CollumnIndex + 1];
						CollumnIndex := CollumnIndex + 2;
					end;
				DoBusyBox(i, n);
				LineIndex := LineIndex + nx4;
			end;
   {    i := 0;   hier keine Symmetrie}
		CollumnIndex := 0;
		for j := 0 to n_1 do begin
				DatVector[j * 2] := FFTDataPtr^[CollumnIndex];
				DatVector[j * 2 + 1] := FFTDataPtr^[CollumnIndex + 1];
				CollumnIndex := CollumnIndex + nx2;
			end;
		FFT(DatVector, n, isign);
		CollumnIndex := 0;
		for j := 0 to n_1 do begin
				FFTDataPtr^[CollumnIndex] := DatVector[j * 2];
				FFTDataPtr^[CollumnIndex + 1] := DatVector[j * 2 + 1];
				CollumnIndex := CollumnIndex + nx2;
			end;
		LineIndex := 2;
		for i := 1 to n_half - 1 do begin
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						DatVector[j * 2] := FFTDataPtr^[LineIndex + CollumnIndex];
						DatVector[j * 2 + 1] := FFTDataPtr^[LineIndex + CollumnIndex + 1];
						CollumnIndex := CollumnIndex + nx2;
					end;
				FFT(DatVector, n, isign);
				CollumnIndex := 0;
				for j := 0 to n_1 - 1 do begin
						FFTDataPtr^[LineIndex + CollumnIndex] := DatVector[j * 2];
						FFTDataPtr^[LineIndex + CollumnIndex + 1] := DatVector[j * 2 + 1];
						FFTDataPtr^[nxnx2 - LineIndex - CollumnIndex] := DatVector[j * 2 + 2];
						FFTDataPtr^[nxnx2 - LineIndex - CollumnIndex + 1] := (-1) * DatVector[j * 2 + 3];
						CollumnIndex := CollumnIndex + nx2;
					end;
				j := n_1;
				FFTDataPtr^[LineIndex + CollumnIndex] := DatVector[j * 2];
				FFTDataPtr^[LineIndex + CollumnIndex + 1] := DatVector[j * 2 + 1];
				FFTDataPtr^[nxnx2 - LineIndex - CollumnIndex] := DatVector[0];
				FFTDataPtr^[nxnx2 - LineIndex - CollumnIndex + 1] := (-1) * DatVector[1];
				DoBusyBox(trunc(n / 2) + i, n);
				LineIndex := LineIndex + 2;
			end;
   {i := n_half;  hier keine Symmetrie;  2*i+n --> LineIndex + CollumnIndex = n + CollumnIndex}
		CollumnIndex := n;
		for j := 0 to n_1 do begin
				DatVector[j * 2] := FFTDataPtr^[CollumnIndex];
				DatVector[j * 2 + 1] := FFTDataPtr^[CollumnIndex + 1];
				CollumnIndex := CollumnIndex + nx2;
			end;
		FFT(DatVector, n, isign);
		CollumnIndex := n;
		for j := 0 to n_1 do begin
				FFTDataPtr^[CollumnIndex] := DatVector[j * 2];
				FFTDataPtr^[CollumnIndex + 1] := DatVector[j * 2 + 1];
				CollumnIndex := CollumnIndex + nx2;
			end;
		HUnlock(handle(FFTArrayH));
		CloseBusyBox;
	end;

	procedure ReFFTHorizontal (ReFFTArrayH: FFTArrayhandle; w, h, isign: longint);
		var
			i, j: longint;
			DatVector: FFTVector;
			LineIndex, CollumnIndex: longint;
			nx2: longint;
			nx4: longint;
			nxnx2: longint;
			n_1: longint;
			n_half: longint;
			fft1, fft2: FFTVector;
			ReFFTDataPtr: FFTArrayPtr;
			DataFile: text;
	begin
		n_1 := w - 1;
		n_half := trunc(h / 2);
		nx2 := w * 2;
		nx4 := w * 4;
		nxnx2 := w * w * 2;
		OpenBusyBox('Processing horizontal Fourier Retrafo');
		LineIndex := 0;
		HLock(Handle(ReFFTArrayH));
		ReFFTDataPtr := ReFFTArrayH^;
		for i := 0 to n_half - 1 do begin
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						DatVector[j * 2] := ReFFTDataPtr^[LineIndex + CollumnIndex] - ReFFTDataPtr^[LineIndex + CollumnIndex + nx2 + 1];
						DatVector[j * 2 + 1] := ReFFTDataPtr^[LineIndex + CollumnIndex + 1] + ReFFTDataPtr^[LineIndex + CollumnIndex + nx2];
						CollumnIndex := CollumnIndex + 2;
					end;
				FFT(DatVector, w, isign);
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						ReFFTDataPtr^[LineIndex + CollumnIndex] := DatVector[j * 2];
						ReFFTDataPtr^[LineIndex + CollumnIndex + 1] := 0.0;
						ReFFTDataPtr^[LineIndex + CollumnIndex + nx2] := DatVector[j * 2 + 1];
						ReFFTDataPtr^[LineIndex + CollumnIndex + nx2 + 1] := 0.0;
						CollumnIndex := CollumnIndex + 2;
					end;
				DoBusyBox(2 * i, h);
				LineIndex := LineIndex + nx4;
			end;
		HUnlock(handle(ReFFTArrayH));
		CloseBusyBox;
	end;



	procedure ReFFT2D (ReFFTArrayH: FFTArrayhandle; n, isign: longint);
		var
			i, j: longint;
			DatVector: FFTVector;
			LineIndex, CollumnIndex: longint;
			nx2: longint;
			nx4: longint;
			nxnx2: longint;
			n_1: longint;
			ReFFTDataPtr: FFTArrayPtr;
			n_half: longint;
			fft1, fft2: FFTVector;
			DataFile: text;
	begin
		n_1 := n - 1;
		n_half := trunc(n / 2);
		nx2 := n * 2;
		nx4 := n * 4;
		nxnx2 := n * n * 2;
		OpenBusyBox('Processing Fourier Retransformation');
    {i := 0; hier keine Symmetrie}
		CollumnIndex := 0;
		HLock(Handle(ReFFTArrayH));
		ReFFTDataPtr := ReFFTArrayH^;
		for j := 0 to n_1 do begin
				DatVector[j * 2] := ReFFTDataPtr^[CollumnIndex];
				DatVector[j * 2 + 1] := ReFFTDataPtr^[CollumnIndex + 1];
				CollumnIndex := CollumnIndex + nx2;
			end;
		FFT(DatVector, n, isign);
		CollumnIndex := 0;
		for j := 0 to n_1 do begin
				ReFFTDataPtr^[CollumnIndex] := DatVector[j * 2];
				ReFFTDataPtr^[CollumnIndex + 1] := DatVector[j * 2 + 1];
				CollumnIndex := CollumnIndex + nx2;
			end;
		LineIndex := 2;
		for i := 1 to n_half - 2 do begin
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						DatVector[j * 2] := ReFFTDataPtr^[LineIndex + CollumnIndex];
						DatVector[j * 2 + 1] := ReFFTDataPtr^[LineIndex + CollumnIndex + 1];
						CollumnIndex := CollumnIndex + nx2;
					end;
				FFT(DatVector, n, isign);
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						ReFFTDataPtr^[LineIndex + CollumnIndex] := DatVector[j * 2];
						ReFFTDataPtr^[LineIndex + CollumnIndex + 1] := DatVector[j * 2 + 1];
						ReFFTDataPtr^[nx2 - LineIndex + CollumnIndex] := DatVector[j * 2];
						ReFFTDataPtr^[nx2 - LineIndex + CollumnIndex + 1] := (-1) * DatVector[j * 2 + 1];
						CollumnIndex := CollumnIndex + nx2;
					end;
				DoBusyBox(i, n);
				LineIndex := LineIndex + 2;
			end;
		for i := n_half - 1 to n_half + 1 do begin
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						DatVector[j * 2] := ReFFTDataPtr^[LineIndex + CollumnIndex];
						DatVector[j * 2 + 1] := ReFFTDataPtr^[LineIndex + CollumnIndex + 1];
						CollumnIndex := CollumnIndex + nx2;
					end;
				FFT(DatVector, n, isign);
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						ReFFTDataPtr^[LineIndex + CollumnIndex] := DatVector[j * 2];
						ReFFTDataPtr^[LineIndex + CollumnIndex + 1] := DatVector[j * 2 + 1];
						CollumnIndex := CollumnIndex + nx2;
					end;
				DoBusyBox(i, n);
				LineIndex := LineIndex + 2;
			end;
		LineIndex := 0;
		for i := 0 to n_half - 1 do begin
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						DatVector[j * 2] := ReFFTDataPtr^[LineIndex + CollumnIndex] - ReFFTDataPtr^[LineIndex + CollumnIndex + nx2 + 1];
						DatVector[j * 2 + 1] := ReFFTDataPtr^[LineIndex + CollumnIndex + 1] + ReFFTDataPtr^[LineIndex + CollumnIndex + nx2];
						CollumnIndex := CollumnIndex + 2;
					end;
				FFT(DatVector, n, isign);
				CollumnIndex := 0;
				for j := 0 to n_1 do begin
						ReFFTDataPtr^[LineIndex + CollumnIndex] := DatVector[j * 2];
						ReFFTDataPtr^[LineIndex + CollumnIndex + 1] := 0.0;
						ReFFTDataPtr^[LineIndex + CollumnIndex + nx2] := DatVector[j * 2 + 1];
						ReFFTDataPtr^[LineIndex + CollumnIndex + nx2 + 1] := 0.0;
						CollumnIndex := CollumnIndex + 2;
					end;
				DoBusyBox(trunc(n / 2) + i, n);
				LineIndex := LineIndex + nx4;
			end;
		HUnlock(handle(ReFFTArrayH));
		CloseBusyBox;
	end;



	procedure ReadPict (FFTArrayH: FFTArrayHandle; n: longint);
		var
			i, j: longint;
			Line: LineType;
			FFTDataPtr: FFTArrayPtr;
			LineIndex, CollumnIndex: longint;
			nx2: longint;
	begin
		nx2 := n * 2;
		LineIndex := 0;
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		for i := 0 to (n - 1) do begin
				GetLine(0, i, n, Line);
				CollumnIndex := 0;
				for j := 0 to (n - 1) do begin
						FFTDataPtr^[LineIndex + CollumnIndex] := Line[j];
						FFTDataPtr^[LineIndex + CollumnIndex + 1] := 0.0;
						CollumnIndex := CollumnIndex + 2;
					end;
				LineIndex := LineIndex + nx2;
			end;
		HUnLock(Handle(FFTArrayH));
	end;


	procedure ReadPictandScale (FFTArrayH: FFTArrayHandle; n: longint);
		var
			i, j, sum: longint;
			Line: LineType;
			FFTDataPtr: FFTArrayPtr;
			LineIndex, CollumnIndex: longint;
			nx2: longint;
			HelpVariable, sum2: real;
	begin
		nx2 := n * 2;
		sum := 0;
		for i := 0 to (n - 1) do begin
				GetLine(0, i, n, Line);
				for j := 0 to (n - 1) do begin
						sum := sum + longint(Line[j]);
					end;
			end;
		sum := sum div longint(n * n);
		LineIndex := 0;
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		for i := 0 to (n - 1) do begin
				GetLine(0, i, n, Line);
				CollumnIndex := 0;
				for j := 0 to (n - 1) do begin
						FFTDataPtr^[LineIndex + CollumnIndex] := Line[j] - sum;
						FFTDataPtr^[LineIndex + CollumnIndex + 1] := 0.0;
						CollumnIndex := CollumnIndex + 2;
					end;
				LineIndex := LineIndex + nx2;
			end;
		LineIndex := 0;
		sum2 := 0;
		for i := 0 to (n - 1) do begin
				CollumnIndex := 0;
				for j := 0 to (n - 1) do begin
						HelpVariable := FFTDataPtr^[LineIndex + CollumnIndex];
						sum2 := sum2 + HelpVariable * HelpVariable;
						CollumnIndex := CollumnIndex + 2;
					end;
				LineIndex := LineIndex + nx2;
			end;
		LineIndex := 0;
		if sum2 = 0 then
			sum2 := 1;
		for i := 0 to (n - 1) do begin
				CollumnIndex := 0;
				for j := 0 to (n - 1) do begin
						FFTDataPtr^[LineIndex + CollumnIndex] := FFTDataPtr^[LineIndex + CollumnIndex] * n / sqrt(sum2);
						CollumnIndex := CollumnIndex + 2;
					end;
				LineIndex := LineIndex + nx2;
			end;


		HUnLock(Handle(FFTArrayH));
	end;


	procedure ReadPictHorizontal (FFTArrayH: FFTArrayHandle; n, height: longint);
		var
			i, j: longint;
			Line: LineType;
			FFTDataPtr: FFTArrayPtr;
			LineIndex, CollumnIndex: longint;
			nx2: longint;
	begin
		nx2 := n * 2;
		LineIndex := 0;
		for i := 0 to (height - 1) do begin
				GetLine(0, i, n, Line);
				CollumnIndex := 0;
				HLock(Handle(FFTArrayH));
				FFTDataPtr := FFTArrayH^;
				for j := 0 to (n - 1) do begin
						FFTDataPtr^[LineIndex + CollumnIndex] := Line[j];
						FFTDataPtr^[LineIndex + CollumnIndex + 1] := 0.0;
						CollumnIndex := CollumnIndex + 2;
					end;
				LineIndex := LineIndex + nx2;
			end;
		HUnLock(Handle(FFTArrayH));
	end;

	procedure ReadPictAndScaleHorizontal (FFTArrayH: FFTArrayHandle; n, height: longint);
		var
			i, j, sum: longint;
			Line, line2: LineType;
			LineIndex, CollumnIndex: longint;
			nx2: longint;
			FFTDataPtr: FFTArrayPtr;
			HelpVariable, sum2: real;
	begin
		nx2 := n * 2;
		LineIndex := 0;
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		for i := 0 to (height - 1) do begin
				GetLine(0, i, n, Line);
				sum := 0;
				for j := 0 to (n - 1) do
					sum := sum + longint(Line[j]);
				sum := sum div n;

				CollumnIndex := 0;
				for j := 0 to (n - 1) do begin
						FFTDataPtr^[LineIndex + CollumnIndex] := Line[j] - sum;
						FFTDataPtr^[LineIndex + CollumnIndex + 1] := 0.0;
						CollumnIndex := CollumnIndex + 2;
					end;
				LineIndex := LineIndex + nx2;
			end;
		LineIndex := 0;
		for i := 0 to (height - 1) do begin
				sum2 := 0;
				CollumnIndex := 0;
				for j := 0 to (n - 1) do begin
						HelpVariable := FFTDataPtr^[LineIndex + CollumnIndex];
						sum2 := sum2 + HelpVariable * HelpVariable;
						CollumnIndex := CollumnIndex + 2;
					end;
				CollumnIndex := 0;
				if sum2 = 0 then
					sum2 := 1;
				for j := 0 to (n - 1) do begin
						FFTDataPtr^[LineIndex + CollumnIndex] := FFTDataPtr^[LineIndex + CollumnIndex] / sqrt(sum2) * sqrt(n);
						CollumnIndex := CollumnIndex + 2;
					end;
				LineIndex := LineIndex + nx2;
			end;
		HUnlock(Handle(FFTArrayH))
	end;



	procedure WriteRealPict (ReFFTArrayH: FFTArrayHandle; w, h: longint);
		var
			i, j: longint;
			LineIndex, CollumnIndex: longint;
			nx2: longint;
			Line: LineType;
			ReFFTDataPtr: FFTArrayPtr;
	begin
		nx2 := w * 2;
		LineIndex := 0;
		HLock(Handle(ReFFTArrayH));
		ReFFTDataPtr := ReFFTArrayH^;
		for i := 0 to (h - 1) do begin
				CollumnIndex := 0;
				for j := 0 to (w - 1) do begin
						Line[j] := trunc(abs(ReFFTDataPtr^[LineIndex + CollumnIndex]));
            { Da Im=0  trunc(  abs (Re)  )  statt  trunc(  sqrt( sqr(Re) + sqr(Im) )  );}
						CollumnIndex := CollumnIndex + 2;
					end;
				PutLine(0, i, w, Line);
				LineIndex := LineIndex + nx2;
			end;
		HUnlock(Handle(ReFFTArrayH))
	end;
	procedure WriteMaxMin (Max, Min: real);
		var
			tPort: GrafPtr;

	begin
		GetPort(tPort);
		SetPort(ResultsWindow);
		TextSize(9);
		TextFont(Monaco);
		TextMode(SrcCopy);
		MoveTo(2, 140);
		DrawBString('Max: ');
		Drawreal(Max, 1, 3);
		DrawBString('  Min: ');
		Drawreal(Min, 1, 3);
		SetPort(tPort);
	end;


	procedure WriteAndScaleRealPict2D (ReFFTArrayH: FFTArrayHandle; w, h: longint; ScalingFactor: real; Button: CorrBut);
		var
			i, j, Temp: longint;
			LineIndex, CollumnIndex: longint;
			nx2: longint;
			Line: LineType;
			max, min: real;
			FFTDataPtr: FFTArrayPtr;
	begin
		HLock(Handle(ReFFTArrayH));
		FFTDataPtr := ReFFTArrayH^;
		nx2 := w * 2;
		max := 0;
		min := 0;
		LineIndex := 0;
		for i := 0 to (h - 1) do begin
				CollumnIndex := 0;
				for j := 0 to (w - 1) do begin
						if FFTDataPtr^[LineIndex + CollumnIndex] > max then
							max := FFTDataPtr^[LineIndex + CollumnIndex];
						if FFTDataPtr^[LineIndex + CollumnIndex] < min then
							min := FFTDataPtr^[LineIndex + CollumnIndex];
						CollumnIndex := CollumnIndex + 2;
					end;
				LineIndex := LineIndex + nx2;
			end;
		WriteMaxMin(max, min);
		if max < -min then
			max := -min;
		if Button = Fixed then
			max := Scalingfactor;
		max := 128 / (max);

		LineIndex := 0;
		for i := 0 to (h div 2 - 1) do begin
				CollumnIndex := 0;
				for j := 0 to (w div 2 - 1) do begin
						Temp := 128 - trunc((FFTDataPtr^[(LineIndex + nx2 * h div 2) + CollumnIndex]) * max);
						if Temp > 255 then
							Temp := 255;
						if Temp < 0 then
							Temp := 0;
						Line[j + w div 2] := Temp;
						CollumnIndex := CollumnIndex + 2;
					end;
				for j := w div 2 to (w - 1) do begin
						Temp := 128 - trunc((FFTDataPtr^[(LineIndex + nx2 * h div 2) + CollumnIndex]) * max);
						if Temp > 255 then
							Temp := 255;
						if Temp < 0 then
							Temp := 0;
						Line[j - w div 2] := Temp;
						CollumnIndex := CollumnIndex + 2;
					end;

				PutLine(0, i, w, Line);
				LineIndex := LineIndex + nx2;
			end;
		LineIndex := nx2 * h div 2;
		for i := h div 2 to (h - 1) do begin
				CollumnIndex := 0;
				for j := 0 to (w div 2 - 1) do begin
						Temp := 128 - trunc((FFTDataPtr^[(LineIndex - nx2 * h div 2) + CollumnIndex]) * max);
						if Temp > 255 then
							Temp := 255;
						if Temp < 0 then
							Temp := 0;
						Line[j + w div 2] := Temp;
						CollumnIndex := CollumnIndex + 2;
					end;
				for j := w div 2 to (w - 1) do begin
						Temp := 128 - trunc((FFTDataPtr^[(LineIndex - nx2 * h div 2) + CollumnIndex]) * max);
						if Temp > 255 then
							Temp := 255;
						if Temp < 0 then
							Temp := 0;
						Line[j - w div 2] := Temp;
						CollumnIndex := CollumnIndex + 2;
					end;

				PutLine(0, i, w, Line);
				LineIndex := LineIndex + nx2;
			end;

		HUnLock(Handle(ReFFTArrayH));
	end;

	procedure WriteAndScaleRealPictHorizontal (ReFFTArrayH: FFTArrayHandle; w, h: longint; ScalingFactor: real; Button: CorrBut);
		var
			i, j, Temp: longint;
			LineIndex, CollumnIndex: longint;
			nx2: longint;
			Line: LineType;
			max, min: real;
			FFTDataPtr: FFTArrayPtr;
	begin
		HLock(Handle(ReFFTArrayH));
		FFTDataPtr := ReFFTArrayH^;
		nx2 := w * 2;
		max := 0;
		min := 0;
		LineIndex := 0;
		for i := 0 to (h - 1) do begin
				CollumnIndex := 0;
				for j := 0 to (w - 1) do begin
						if FFTDataPtr^[LineIndex + CollumnIndex] > max then
							max := FFTDataPtr^[LineIndex + CollumnIndex];
						if FFTDataPtr^[LineIndex + CollumnIndex] < min then
							min := FFTDataPtr^[LineIndex + CollumnIndex];
						CollumnIndex := CollumnIndex + 2;
					end;
				LineIndex := LineIndex + nx2;
			end;
		WriteMaxMin(max, min);
		if max < -min then
			max := -min;
		if Button = Fixed then
			max := Scalingfactor;
		max := 128 / (max);
		LineIndex := 0;
		for i := 0 to (h - 1) do begin
				CollumnIndex := 0;
				for j := 0 to (w div 2 - 1) do begin
						Temp := 128 - trunc((FFTDataPtr^[LineIndex + CollumnIndex]) * max);
						if Temp > 255 then
							Temp := 255;
						if Temp < 0 then
							Temp := 0;
						Line[j + w div 2] := Temp;
						CollumnIndex := CollumnIndex + 2;
					end;
				for j := w div 2 to (w - 1) do begin
						Temp := 128 - trunc((FFTDataPtr^[LineIndex + CollumnIndex]) * max);
						if Temp > 255 then
							Temp := 255;
						if Temp < 0 then
							Temp := 0;
						Line[j - w div 2] := Temp;
						CollumnIndex := CollumnIndex + 2;
					end;

				PutLine(0, i, w, Line);
				LineIndex := LineIndex + nx2;
			end;
		HUnLock(Handle(ReFFTArrayH));
	end;


	function FindMax (FFTArrayH: FFTArrayHandle; w, h: longint): real;
		var
			i: longint;
			max: real;
			FFTDataPtr: FFTArrayPtr;
	begin
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		max := abs(FFTDataPtr^[3]);
{ skip the first entry which is the dc part }
		for i := 3 to longint(2 * w - 1) * longint(h) do begin
				if max <= abs(FFTDataPtr^[i]) then
					max := abs(FFTDataPtr^[i]);
			end;
		FindMax := max;
		HUnLock(Handle(FFTArrayH));
	end;

	function FindMin (FFTArrayH: FFTArrayHandle; w, h: longint): real;
		var
			i: longint;
			min, auxnum: real;
			FFTDataPtr: FFTArrayPtr;
	begin
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		min := 1;
		for i := 0 to longint(2 * w - 1) * longint(h) do begin
				auxnum := abs(FFTDataPtr^[i]);
				if (auxnum > 0.0) and (min > auxnum) then
					min := auxnum;
			end;
		FindMin := min;
		HUnLock(Handle(FFTArrayH));
	end;

	function FindHorizontalMax (FFTArrayH: FFTArrayHandle; n, height: longint): real;
		var
			line, column, offset: longint;
			FFTDataPtr: FFTArrayPtr;
			max: real;
	begin
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		max := abs(FFTArrayH^^[3]);
		offset := 0;
{ skip the first entry which is the dc part }
		for line := 0 to height - 1 do begin
				for column := 3 to 2 * n - 1 do begin
						if max <= abs(FFTDataPtr^[column + offset]) then
							max := abs(FFTDataPtr^[column + offset]);
					end;
				offset := offset + 2 * n;
			end;
		FindHorizontalMax := max;
		HUnLock(Handle(FFTArrayH));
	end;

	function FFT2Color (DataPtr: FFTArrayPtr; a, b: real; index: longint): integer;
		var
			aux: real;
			num: integer;
	begin
		aux := 180 * (ln(sqr(DataPtr^[(index) * 2]) + sqr(DataPtr^[(index) * 2 + 1])) - b) / (a - b);
{ nur 128 *sqrt(2) , da sqrt nicht gemacht }
		num := trunc(aux);
{ scale a liitle bit }
		if num > 255 then
			num := 255
		else if num < 0 then
			num := 0;
		FFT2Color := num;
	end;

	procedure WriteFFTPict (FFTArrayH: FFTArrayHandle; n: longint);
		var
			i, j, LineOffset: longint;
			aNumber: integer;
			Line, Line2: LineType;
			n_Half: longint;
			a, b: real;
			FFTDataPtr: FFTArrayPtr;
	begin
		a := 1;
		if a <> 0 then
			a := ln(FindMax(FFTArrayH, n, N));
		b := ln(FindMin(FFTArrayH, n, N));
		n_half := trunc(n / 2);
		LineOffset := 0;
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;

{ the first line }
		for j := 0 to n_half - 1 do
			Line[n_half + j] := FFT2Color(FFTDataPtr, a, b, j);
		for j := n_half to (n - 1) do
			Line[j - n_half] := FFT2Color(FFTDataPtr, a, b, j);
		LineOffset := LineOffset + n;
		PutLine(0, n_half, n, Line);

		for i := 1 to n_half - 1 do begin
				for j := 0 to n_half - 1 do
					Line[n_half + j] := FFT2Color(FFTDataPtr, a, b, LineOffset + j);

				for j := n_half to (n - 1) do
					Line[j - n_half] := FFT2Color(FFTDataPtr, a, b, LineOffset + j);
				PutLine(0, n_half + i, n, Line);
				for j := 1 to n - 1 do
					Line2[n - j] := Line[j];
				Line2[0] := Line[0];
				LineOffset := LineOffset + n;
				PutLine(0, n_half - i, n, Line2);
			end;

{ there is one line left }
		for j := 0 to n_half - 1 do
			Line[n_half + j] := FFT2Color(FFTDataPtr, a, b, LineOffset + j);
		for j := n_half to (n - 1) do
			Line[j - n_half] := FFT2Color(FFTDataPtr, a, b, LineOffset + j);
		PutLine(0, 0, n, Line);
		HUnLock(Handle(FFTArrayH));
		UpDatePicWindow;
	end;

	procedure WriteHorizontalFFTPict (FFTArrayH: FFTArrayHandle; n, height: longint);
		var
			i, j, LineOffset: longint;
			aNumber: integer;
			Line: LineType;
			n_Half: longint;
			a, b: real;
			FFTDataPtr: FFTArrayPtr;
	begin
		a := 1;
		if a <> 0 then
			a := ln(FindMax(FFTArrayH, n, height));
		b := ln(FindMin(FFTArrayH, n, height));
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;

		n_half := trunc(n / 2);
		LineOffset := 0;
		for i := 0 to height - 1 do begin
				for j := 0 to n_half - 1 do
					Line[n_half + j] := FFT2Color(FFTDataPtr, a, b, LineOffset + j);

				for j := n_half to (n - 1) do
					Line[j - n_half] := FFT2Color(FFTDataPtr, a, b, LineOffset + j);
				PutLine(0, i, n, Line);
				LineOffset := LineOffset + n;
			end;
		HUnLock(Handle(FFTArrayH));
	end;

	procedure WriteHorizontalFFTPhases (FFTArrayH: FFTArrayHandle; n, height: longint);
		var
			i, j: longint;
			aNumber: integer;
			Line: LineType;
			n_Half: longint;
			FFTDataPtr: FFTArrayPtr;
			help: real;
	begin
		n_half := trunc(n / 2);
		HLock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		for i := 0 to height - 1 do begin
				for j := 0 to n_half do begin
						aNumber := trunc((0.5 * Pi + arctan(FFTDataPtr^[(i * n + j) * 2 + 1] / FFTDataPtr^[(i * n + j) * 2])) * 256 / PI);
						if aNumber > 255 then
							Line[j] := 255
						else
							Line[j] := Band(aNumber, 255);
					end;
				PutLine(n_half, i, n_half, Line);
				for j := n_half to (n - 1) do begin
						aNumber := trunc((0.5 * Pi + arctan(FFTDataPtr^[(i * n + j) * 2 + 1] / FFTDataPtr^[(i * n + j) * 2])) * 256 / PI);
						if aNumber > 255 then
							Line[j - n_half] := 255
						else
							Line[j - n_half] := Band(aNumber, 255);
					end;
				PutLine(0, i, n_half, Line);
			end;
		HUnLock(Handle(FFTArrayH));
	end;


	procedure SaveFFTPicture (SavedPictH: SavedPictHandle; n: longint);
		var
			i, j: longint;
			Line: LineType;
			LineIndex: longint;
			SavedDataptr: SavedPictPtr;
	begin
		LineIndex := 0;
		Hlock(Handle(SavedPictH));
		SavedDataPtr := SavedPictH^;
		for i := 0 to (n - 1) do begin
				GetLine(0, i, n, Line);
				for j := 0 to (n - 1) do begin
						SavedDataPtr^[LineIndex + j] := Line[j];
					end;
				LineIndex := LineIndex + n;
			end;
		HUnlock(Handle(SavedPictH));
	end;



	procedure RestoreFFTPicture (SavedPictH: SavedPictHandle; n: longint);
		var
			i, j: longint;
			LineIndex: longint;
			Line: LineType;
			SaveInfo: InfoPtr;
			SavedDataptr: SavedPictPtr;
	begin
		SaveInfo := Info;
		with FFTInfo do
			Info := FFTPict_Info;
		LineIndex := 0;
		Hlock(Handle(SavedPictH));
		SavedDataPtr := SavedPictH^;
		for i := 0 to (n - 1) do begin
				for j := 0 to (n - 1) do begin
						Line[j] := SavedDataptr^[LineIndex + j];
					end;
				PutLine(0, i, n, Line);
				LineIndex := LineIndex + n;
			end;
		UpdatePicWindow;
		Info := SaveInfo;
		HUnlock(Handle(SavedPictH));
	end;



	procedure CopyFFTArray (FFTArrayH, ReFFTArrayH: FFTArrayHandle; w, h: longint);
		var
			size: longint;
	begin
		Hlock(handle(FFTArrayH));
		Hlock(handle(ReFFTArrayH));
		size := 2 * w * h * sizeof(real);
		BlockMove(ptr(FFTArrayH^), ptr(ReFFTArrayH^), size);
		HUnlock(handle(FFTArrayH));
		HUnlock(handle(ReFFTArrayH));
	end;

	procedure Do2DFFT;
		var
			DatVector: FFTVector;
			w, h: longint;
	begin
		if DuplicateFFT(false, w, h) then begin
				with info^ do begin
						ShowWatch;
						ReadPict(FFTArrayH, w);
						FFT2D(FFTArrayH, w, 1);
						WriteFFTPict(FFTArrayH, w);
						SaveFFTPicture(SavedPictH, w);
						with FFTInfo do begin
								FFT_Done := true;
								FFTPict_Info := Info;
							end;
					end;
			end
		else
			;
	end;

	procedure DoHorizontalFFT;
		var
			DatVector: FFTVector;
			w, h: longint;
	begin
		if not info^.RoiShowing then
			SelectAll(false);
		with info^.RoiRect do begin
				if (longint(right - left) * longint(bottom - top)) > (NNNN div 2) then begin
						putmessage('selection must be smaller !');
						exit(DoHorizontalFFT);
					end;
			end;
		if DuplicateHorizontalFFT(false, w, h) then begin
				with info^ do begin
						ShowWatch;
						ReadPictHorizontal(FFTArrayH, w, h);
						FFTHorizontal(FFTArrayH, w, h, 1);
						WriteHorizontalFFTPict(FFTArrayH, w, h);
						with FFTInfo do begin
								FFT_Done := true;
								FFTPict_Info := Info;
							end;
					end;
			end
		else
			;
	end;


	procedure GetValue (FFTDataPtr: FFTArrayPtr; n, i, j: longint; var Re, Im: real);
		var
			n_half: longint;
	begin
		n_half := trunc(n / 2);
		if j < n_half then begin
				if i < n_half then begin
						Re := FFTDataptr^[((i + n_half) * n + (j + n_half)) * 2];
						Im := FFTDataptr^[((i + n_half) * n + (j + n_half)) * 2 + 1];
					end
				else begin
						Re := FFTDataptr^[((i - n_half) * n + (j + n_half)) * 2];
						Im := FFTDataptr^[((i - n_half) * n + (j + n_half)) * 2 + 1];
					end;
			end
		else begin
				if i < n_half then begin
						Re := FFTDataptr^[((i + n_half) * n + (j - n_half)) * 2];
						Im := FFTDataptr^[((i + n_half) * n + (j - n_half)) * 2 + 1];
					end
				else begin
						Re := FFTDataptr^[((i - n_half) * n + (j - n_half)) * 2];
						Im := FFTDataptr^[((i - n_half) * n + (j - n_half)) * 2 + 1];
					end;
			end;
	end;


	procedure WriteValue (ReFFTDataPtr: FFTArrayPtr; n, i, j: longint; Re, Im: real);
		var
			n_half: longint;
	begin
		n_half := trunc(n / 2);
		if j < n_half then begin
				if i < n_half then begin
						ReFFTDataptr^[((i + n_half) * n + (j + n_half)) * 2] := Re;
						ReFFTDataptr^[((i + n_half) * n + (j + n_half)) * 2 + 1] := Im;
					end
				else begin
						ReFFTDataptr^[((i - n_half) * n + (j + n_half)) * 2] := Re;
						ReFFTDataptr^[((i - n_half) * n + (j + n_half)) * 2 + 1] := Im;
					end;
			end
		else begin
				if i < n_half then begin
						ReFFTDataptr^[((i + n_half) * n + (j - n_half)) * 2] := Re;
						ReFFTDataptr^[((i + n_half) * n + (j - n_half)) * 2 + 1] := Im;
					end
				else begin
						ReFFTDataptr^[((i - n_half) * n + (j - n_half)) * 2] := Re;
						ReFFTDataptr^[((i - n_half) * n + (j - n_half)) * 2 + 1] := Im;
					end;
			end;
	end;



	procedure WriteChanges (SavedPictH: SavedPictHandle; ReFFTArrayH: FFTArrayHandle; n: longint);
		var
			i, j: longint;
			LineIndex: longint;
			Line: LineType;
			n_half: longint;
			ReFFTDataPtr: FFTArrayPtr;
			SavedDataPtr: SavedPictPtr;
			SaveInfo: InfoPtr;
	begin
		n_half := trunc(n / 2);
		SaveInfo := Info;
		with FFTInfo do
			Info := FFTPict_Info;
		GetLine(0, 0, n, Line);
		Hlock(Handle(ReFFTArrayH));
		ReFFTDataPtr := ReFFTArrayH^;
		Hlock(Handle(SavedPictH));
		SavedDataPtr := SavedPictH^;
		for j := 0 to (n - 1) do begin
				if Line[j] <> SavedDataPtr^[j] then
					WriteValue(ReFFTDataPtr, n, 0, j, 0.0, 0.0);
			end;
		GetColumn(0, 0, n, Line);
		for i := 0 to (n - 1) do begin
				if Line[i] <> SavedDataPtr^[i * n] then
					WriteValue(ReFFTDataPtr, n, i, 0, 0.0, 0.0);
			end;
		for i := 1 to (n_half - 1) do begin
				GetLine(0, i, n, Line);
				if Line[0] <> SavedDataPtr^[i * n] then
					WriteValue(ReFFTDataPtr, n, i, 0, 0.0, 0.0);
				for j := 1 to (n - 1) do begin
						if Line[j] <> SavedDataPtr^[i * n + j] then begin
								WriteValue(ReFFTDataPtr, n, i, j, 0.0, 0.0);
								WriteValue(ReFFTDataPtr, n, n - i, n - j, 0.0, 0.0);
							end;
					end;
			end;
		GetLine(0, n_half, n, Line);
		for j := 0 to (n - 1) do begin
				if Line[j] <> SavedDataPtr^[n_half * n + j] then
					WriteValue(ReFFTDataPtr, n, n_half, j, 0.0, 0.0);
			end;
		Info := SaveInfo;
		HUnlock(Handle(ReFFTArrayH));
		HUnlock(Handle(SavedPictH));
	end;



	procedure SmoothPicture (FFTArrayH, ReFFTArrayH: FFTArrayHandle; n: longint);
		var
			i, j: longint;
			FFTDataPtr, ReFFTDataPtr: FFTArrayPtr;
			Line: LineType;
			SaveInfo: InfoPtr;
			Re, Im: real;
			n_half: longint;
	begin
		n_half := trunc(n / 2);
		SaveInfo := Info;
		with FFTInfo do
			Info := FFTPict_Info;
		GetLine(0, 0, n, Line);
		Hlock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		Hlock(Handle(ReFFTArrayH));
		ReFFTDataPtr := ReFFTArrayH^;
		for j := 0 to (n - 2) do begin
				if ((Line[j] = 0) and (Line[j + 1] <> 0)) then begin
						GetValue(FFTDataPtr, n, 0, j + 1, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, 0, j, Re, Im);
					end;
				if ((Line[n - j - 1] = 0) and (Line[n - j - 2] <> 0)) then begin
						GetValue(FFTDataPtr, n, 0, n - j - 2, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, 0, n - j - 1, Re, Im);
					end;
			end;
		for i := 1 to (n_half - 1) do begin
				GetLine(0, i, n, Line);
				if ((Line[0] = 0) and (Line[1] <> 0)) then begin
						GetValue(FFTDataPtr, n, i, 1, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, i, 0, Re, Im);
					end;
				if ((Line[n - 1] = 0) and (Line[n - 2] <> 0)) then begin
						GetValue(FFTDataPtr, n, i, n - 2, Re, Im);
						WriteValue(ReFFTDataPtr, n, i, n - 1, Re, Im);
						WriteValue(ReFFTDataPtr, n, n - i, 1, Re, Im);
					end;
				for j := 1 to (n - 2) do begin
						if ((Line[j] = 0) and (Line[j + 1] <> 0)) then begin
								GetValue(FFTDataPtr, n, i, j + 1, Re, Im);
								Re := 0.5 * Re;
								Im := 0.5 * Im;
								WriteValue(ReFFTDataPtr, n, i, j, Re, Im);
								WriteValue(ReFFTDataPtr, n, n - i, n - j, Re, Im);
							end;
						if ((Line[n - j - 1] = 0) and (Line[n - j - 2] <> 0)) then begin
								GetValue(FFTDataPtr, n, i, n - j - 2, Re, Im);
								Re := 0.5 * Re;
								Im := 0.5 * Im;
								WriteValue(ReFFTDataPtr, n, i, n - j - 1, Re, Im);
								WriteValue(ReFFTDataPtr, n, n - i, j + 1, Re, Im);
							end;
					end;
			end;
		GetLine(0, n_half, n, Line);
		for j := 0 to (n - 2) do begin
				if ((Line[j] = 0) and (Line[j + 1] <> 0)) then begin
						GetValue(FFTDataPtr, n, n_half, j + 1, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, n_half, j, Re, Im);
					end;
				if ((Line[n - j - 1] = 0) and (Line[n - j - 2] <> 0)) then begin
						GetValue(FFTDataPtr, n, n_half, n - j - 2, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, n_half, n - j - 1, Re, Im);
					end;
			end;
		GetColumn(0, 0, n, Line);
		for j := 0 to (n - 2) do begin
				if ((Line[j] = 0) and (Line[j + 1] <> 0)) then begin
						GetValue(FFTDataPtr, n, j + 1, 0, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, j, 0, Re, Im);
					end;
				if ((Line[n - j - 1] = 0) and (Line[n - j - 2] <> 0)) then begin
						GetValue(FFTDataPtr, n, n - j - 2, 0, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, n - j - 1, 0, Re, Im);
					end;
			end;
		for i := 1 to (n_half - 1) do begin
				GetColumn(i, 0, n, Line);
				if ((Line[0] = 0) and (Line[1] <> 0)) then begin
						GetValue(FFTDataPtr, n, 1, i, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, 0, i, Re, Im);
					end;
				if ((Line[n - 1] = 0) and (Line[n - 2] <> 0)) then begin
						GetValue(FFTDataPtr, n, n - 2, i, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, n - 1, i, Re, Im);
						WriteValue(ReFFTDataPtr, n, 1, n - i, Re, Im);
					end;
				for j := 1 to (n - 2) do begin
						if ((Line[j] = 0) and (Line[j + 1] <> 0)) then begin
								GetValue(FFTDataPtr, n, i, j + 1, Re, Im);
								Re := 0.5 * Re;
								Im := 0.5 * Im;
								WriteValue(ReFFTDataPtr, n, j, i, Re, Im);
								WriteValue(FFTDataPtr, n, n - j, n - i, Re, Im);
							end;
						if ((Line[n - j - 1] = 0) and (Line[n - j - 2] <> 0)) then begin
								GetValue(FFTDataPtr, n, n - j - 2, i, Re, Im);
								Re := 0.5 * Re;
								Im := 0.5 * Im;
								WriteValue(ReFFTDataPtr, n, n - j - 1, i, Re, Im);
								WriteValue(ReFFTDataPtr, n, j + 1, n - i, Re, Im);
							end;
					end;
			end;
		GetColumn(n_half, 0, n, Line);
		for j := 0 to (n - 2) do begin
				if ((Line[j] = 0) and (Line[j + 1] <> 0)) then begin
						GetValue(FFTDataPtr, n, j + 1, n_half, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, j, n_half, Re, Im);
					end;
				if ((Line[n - j - 1] = 0) and (Line[n - j - 2] <> 0)) then begin
						GetValue(FFTDataPtr, n, n - j - 2, n_half, Re, Im);
						Re := 0.5 * Re;
						Im := 0.5 * Im;
						WriteValue(ReFFTDataPtr, n, n - j - 1, n_half, Re, Im);
					end;
			end;
		Info := SaveInfo;
		HUnlock(Handle(FFTArrayH));
		HUnlock(Handle(ReFFTArrayH));
	end;



	procedure Do2DBackFFT;
		var
			w, h: longint;
	begin
		ShowWatch;
		with FFTInfo.FFTPict_Info^ do begin
				CopyFFTArray(FFTArrayH, ReFFTArrayH, PixelsPerLine, nlines);
				WriteChanges(SavedPictH, ReFFTArrayH, PixelsPerLine);
				SmoothPicture(FFTArrayH, ReFFTArrayH, PixelsPerLine);
			end;
		if DuplicateFFT(false, w, h) then begin
				ReFFT2D(ReFFTArrayH, w, -1);
				WriteRealPict(ReFFTArrayH, w, h);
			end;
	end;



	procedure Restore;
	begin
		ShowWatch;
		with FFTInfo.FFTPict_Info^ do begin
				CopyFFTArray(FFTArrayH, ReFFTArrayH, PixelsPerLine, nlines);
				RestoreFFTPicture(SavedPictH, PixelsPerLine);
			end;
	end;



	procedure OpenBusyBox;
		var
			Itemtype: integer;
			Item: handle;
			tport: GrafPtr;
	begin
		BusyDialog := GetNewDialog(BusyDLOG, nil, pointer(-1));
		GetDitem(BusyDialog, BusyBox, ItemType, item, BusyRect);
		SetDString(BusyDialog, 1, text);
		DrawDialog(BusyDialog);
		GetPort(tPort);
		SetPort(BusyDialog);
		FrameRect(BusyRect);
		SetPort(tPort);
	end;



	procedure DoBusyBox;
		var
			tport: GrafPtr;
	begin
		if a > 0 then begin
				GetPort(tPort);
				SetPort(BusyDialog);
				BusyRect.right := BusyRect.left + (longint(400) * n div a);
				FillRect(BusyRect, gray);
				SetPort(tPort);
			end;
	end;


	procedure CloseBusyBox;
	begin
		DisposDialog(BusyDialog);
	end;



	procedure DoFlaten (components, n: integer);
		var
			radius, i, j, n_half: longint;
			SaveInfo: InfoPtr;
			Line: LineType;
	begin
		SaveInfo := Info;
		Info := FFTInfo.FFTPict_Info;
		n_half := n div 2;
		with Info^ do begin
				for i := n_half - components to n_half + components do begin
						GetLine(0, i, n, Line);
						for j := n_half - components to n_half + components do begin
								if (i <> n_half) or (j <> n_half) then
									Line[j] := 0;
							end;
						PutLine(0, i, n, Line);
					end;
			end;
		UpdatePicWindow;
		Info := SaveInfo;
	end;

	procedure CutOffNoise (noise, n: integer);
		var
			radius: longint;
			i, j: longint;
			aux, n_half: integer;
			SaveInfo: InfoPtr;
			Line: LineType;
	begin
		SaveInfo := Info;
		Info := FFTInfo.FFTPict_Info;
		n_half := n div 2;
		UpdatePicWindow;
		with Info^ do begin
				radius := trunc((1 - noise / 100) * n / 2);
				for i := 0 to (n - 1) do begin
						if i < n_half then begin
								if i < (n_half - radius) then
									aux := n_half
								else
									aux := n_half - trunc(sqrt(radius * radius - (n_half - i) * (n_half - i)));
							end;
						if i > n_half then begin
								if i > (n_half + radius) then
									aux := n_half
								else
									aux := n_half - trunc(sqrt(radius * radius - (i - n_half) * (i - n_half)));
							end;
						GetLine(0, i, n, Line);
						for j := 0 to aux do
							Line[j] := 0;
						for j := n - aux to n do
							Line[j] := 0;
						PutLine(0, i, n, Line);
					end;
			end;
		UpdatePicWindow;
	end;



	procedure DoFFTFilter;
		var
			mylog: DialogPtr;
			item, i: integer;
	begin
		mylog := GetNewDialog(FilterDLOG, nil, pointer(-1));
		OutlineButton(mylog, FilterOK, 16);
		SetDnum(mylog, NoisePercent, FFTFilterNoise);
		SetDnum(mylog, FlatenText, FFTFilterComponents);

		repeat
			begin
				ModalDialog(nil, item);
				case item of
					NoiseFastUp: 
						begin
							if FFTFilterNoise < 90 then
								FFTFilterNoise := FFTFilterNoise + 10;
							SetDnum(mylog, NoisePercent, FFTFilterNoise);
						end;
					NoiseFastDown: 
						begin
							if FFTFilterNoise >= 10 then
								FFTFilterNoise := FFTFilterNoise - 10;
							SetDnum(mylog, NoisePercent, FFTFilterNoise);
						end;
					NoiseUp: 
						begin
							if FFTFilterNoise < 100 then
								FFTFilterNoise := FFTFilterNoise + 1;
							SetDnum(mylog, NoisePercent, FFTFilterNoise);
						end;
					NoiseDown: 
						begin
							if FFTFilterNoise > 0 then
								FFTFilterNoise := FFTFilterNoise - 1;
							SetDnum(mylog, NoisePercent, FFTFilterNoise);
						end;
					FlatenUp: 
						begin
							if FFTFilterComponents < 50 then
								FFTFilterComponents := FFTFilterComponents + 1;
							SetDnum(mylog, FlatenText, FFTFilterComponents);
						end;
					FlatenDown: 
						begin
							if FFTFilterComponents >= 1 then
								FFTFilterComponents := FFTFilterComponents - 1;
							SetDnum(mylog, FlatenText, FFTFilterComponents);
						end;
					FlatenFastUp: 
						begin
							if FFTFilterComponents < 50 then
								FFTFilterComponents := FFTFilterComponents + 10;
							SetDnum(mylog, FlatenText, FFTFilterComponents);
						end;
					FlatenFastDown: 
						begin
							if FFTFilterComponents >= 10 then
								FFTFilterComponents := FFTFilterComponents - 10;
							SetDnum(mylog, FlatenText, FFTFilterComponents);
						end;
					otherwise
				end;
			end;
		until (item = FilterOk) or (item = FilterCancel);

		if item = FilterOk then begin
				DoFlaten(FFTFilterComponents, n);
				CutOffNoise(FFTFilterNoise, n);
			end;
		DisposDialog(mylog);
	end;

	var
		k_Vector: integer;

	function GetkToSave: integer;
		const
			SaveKDLOG = 2301;
			SaveKOk = 1;
			SaveKCancel = 2;
			SaveKText = 5;
		var
			mylog: DialogPtr;
			item, i: integer;
	begin
		mylog := GetNewDialog(SaveKDLog, nil, pointer(-1));
		OutlineButton(mylog, SaveKOK, 16);
		SetDNum(mylog, SaveKText, k_Vector);
		SelIText(mylog, SaveKText, 0, 32767);
		repeat
			ModalDialog(nil, item);
		until (item = saveKOk) or (item = saveKCancel);
		if item = saveKOk then
			GetkToSave := Round(GetDReal(mylog, SaveKText))
		else
			GetkToSave := -1;
		DisposDialog(mylog);
	end;

	procedure saveFFTof_t (FFTArrayH: FFTArrayHandle; n, height: longint);
		const
			CustomDialogID = 60;
		var
			i, j: longint;
			result: integer;
			Line: LineType;
			n_Half, count: longint;
			help: real;
			Filename, name: Str255;
			ReFFT, ImFFT: array[0..500] of real;
			RefNum, vRefNum: integer;
			buff: Str255;
			where: Point;
			reply: SFReply;
			isSelection: boolean;
			kind: integer;
			FFTDataPtr: FFTArrayPtr;
	begin
		n_half := trunc(n / 2);
		if height - 1 > 499 then begin
				putMessage('arrays in procedure  saveFFTof_t are declared too small');
				exit(saveFFTof_t);
			end;
		k_Vector := GetkToSave;
		if k_Vector > n_half then begin
				putMessage('k must not be greater than  half of the width of FFT-Picture');
				exit(saveFFTof_t);
			end;
		if k_Vector < 0 then begin
				exit(saveFFTof_t);
			end;
		Hlock(Handle(FFTArrayH));
		FFTDataPtr := FFTArrayH^;
		for i := 0 to height - 1 do begin
				j := k_Vector;
				ImFFT[i] := FFTDataPtr^[(i * n + j) * 2 + 1];
				ReFFT[i] := FFTDataPtr^[(i * n + j) * 2];
			end;
		HUnlock(Handle(FFTArrayH));
		with Info^ do {FFTInfo.FFTPict_}
			begin
				GetWTitle(wptr, name);
				FileName := concat('FFT(', StringOf(k_Vector : 0), ')', copy(name, 7, length(name)));
				result := Create(Filename, vRef, '    ', 'TEXT');
				if result <> 0 then begin
						where.v := 50;
						where.h := 50;
						SFPPutFile(Where, Concat('ResC=', StringOf(result : 0)), Filename, nil, reply, CustomDialogID, nil);
						if not reply.good then begin
								exit(saveFFTof_t);
							end;
						vref := reply.vRefNum;
						FileName := reply.fName;
						result := FSDelete(FileName, vref);
						result := Create(Filename, vRef, '    ', 'TEXT');
					end;
				result := FSOpen(FileName, vRef, RefNum);
				buff := concat('IGOR', chr(13), 'WAVES', chr(9), 'Re', copy(name, 7, 8), '_', StringOf(k_Vector : 0), chr(9), 'Im', copy(name, 7, 8), '_', StringOf(k_Vector : 0), chr(13), 'BEGIN', chr(13));
				count := Length(buff);
				result := FSWrite(RefNum, count, Ptr(longint(@buff) + 1));
				for i := 0 to height - 1 do begin
						buff := StringOf(ReFFT[i] : 10 : 5, chr(9), ImFFT[i] : 10 : 5, chr(13));
						count := Length(buff);
						result := FSWrite(RefNum, count, Ptr(longint(@buff) + 1));
					end;
				buff := concat('END', chr(13));
				count := Length(buff);
				result := FSWrite(RefNum, count, Ptr(longint(@buff) + 1));
				result := FSClose(RefNum);
			end;
	end;


	procedure ShowFFTPhases;
		var
			name: str255;
			width, height, hstart, vstart: integer;
			i: integer;
			SaveInfo: InfoPtr;
	begin
		with FFTInfo.FFTPict_Info^ do begin
				width := PixelsPerLine;
				height := nLines;
				GetWTitle(wptr, name);
			end;
		name := concat('Phases of ', name);
		SaveInfo := Info;
{    if NewPicWindow(name, width, height) then}
		saveFFTof_t(FFTArrayH, width, height);
{WriteHorizontalFFTPhases(width, height);}
	end;



end.