unit cppDCT;

//http://www.codeproject.com/Articles/151043/Iterative-Fast-1D-Forvard-DCT

Interface

procedure SlowDCT(var block:array of double; N:integer);
procedure FastDCT(var block:array of double;N:integer);

Implementation

procedure SlowDCT(var block:array of double; N:integer);
  var
	i, j : integer;
  s, PI_by_N,pt:extended;
  tmp:array of double;
  begin
	setLength(tmp,N);
  //block[0]:=1;block[1]:=2;block[2]:=3;block[3]:=4;   //test
  pt:=sqrt(2)/2;
	PI_by_N:=PI/N;
	for i:=0 to N-1 do
	begin
		s:=0;
		for j:=0 to N-1 do s:=s+ block[j]*cos(PI_by_N*i*(j+0.5));
		tmp[i]:=s;
	end;
	block[0]:=tmp[0];
	for i:=1 to N-1 do 	block[i]:=tmp[i]*pt;
	setLength(tmp,0);
end;

procedure FastDCT(var block:array of double;N:integer);
 var
          x,i,j,half_N, prev_half_N, shif_Block, step, shift : integer;
 dN,current_PI_By_N,PI_half,current_PI_half_By_N, step_Phase : double;           //*Front, *Back, *temp; Pointers on
                             Front,Back,temp,Arr,Current_Cos : array of double; //Front, Back, and temp - massives of iteration.
  begin                                                                    //Arr - a massive with doubled size
  setLength(Arr,N shl 1); setLength(Current_Cos,N shr 1);                 //Current_Cos store for precomputing
	Front:=Arr;                                                            // Front pointed on begin of the massive
	step:=0;
	Back:=addr(Arr[N]);//addr(Arr[N]);//Front+N;                                                       // Back pointed on center of the massive
	dN:=N;
	half_N:=N;
	j:=N;
	step:=1;
	prev_half_N:=N;
	for i:=0 to N-1 do 	Front[i]:=block[i];

	while j>1 do                                                 // Cycle of iterations Input Butterfly
	Begin
		half_N:=half_N shr 1;
		current_PI_half_By_N:=PI_half/prev_half_N;
		current_PI_By_N:=0;
		step_Phase:=current_PI_half_By_N*2.0;
		for i:=0 to half_N-1 do                              //Precompute Cosine's coefficients
		begin
			Current_Cos[i]:=0.5/cos(current_PI_By_N+current_PI_half_By_N);
			current_PI_By_N :=current_PI_By_N+step_Phase;
		end;{for i}
		shif_Block:=0;
		for x:=0 to step-1 do
		begin
			for i:=0 to half_N-1 do
			begin
				shift:=shif_Block+prev_half_N-i-1;
				Back[shif_Block+i]:=Front[shif_Block+i]+Front[shift];
				Back[shif_Block+half_N+i]:=(Front[shif_Block+i]-Front[shift])*Current_Cos[i];
			end;{for i}
			shif_Block:= shif_Block+prev_half_N;
		end;{for x}
		temp:=Front;
		Front:=Back;
		Back:=temp;
		j:=j shr 1;
		step:=step shl 1;
		prev_half_N:=half_N;
	End; {while j>1}

	step:=step shr 2;
	j:=j shl 1;
	prev_half_N:=prev_half_N shl 1;
	half_N:=prev_half_N;
	prev_half_N:=prev_half_N shl 1;

	while(j<N) do                                                            // Cycle of Out ButterFly
	Begin
		shif_Block:=0;
		for x:=0 to step-1 do
		begin
			for i:=0 to half_N-2 do
			begin
				Back[shif_Block+(i shl 1)]:=Front[shif_Block+i];
				Back[shif_Block+(i shl 1)+1]:=Front[shif_Block+half_N+i]+Front[shif_Block+half_N+i+1];
			end;{for i}
			Back[shif_Block+((half_N-1) shl 1)]:=Front[shif_Block+half_N-1];
			Back[shif_Block+(half_N shl 1)-1]:=Front[shif_Block+(half_N shl 1)-1];
			shif_Block:=shif_Block+prev_half_N;
		end;{for x}
		temp:=Front;
		Front:=Back;
		Back:=temp;
		j:=j shl 1;
		step:=step shr 1;
		half_N:=prev_half_N;
		prev_half_N:=prev_half_N shl 1;
	End; {while j<N}
	for i:=0 to N-1 do block[i]:=Front[i];  // Unload DCT coefficients
	block[0]:=block[0]/dN;                  // Compute DC.
	setLength(Current_Cos,0);
	setLength(Arr,0);
end;
end.
//void DCT::FastDCT(double *block, int N)
//begin
//	int i, j, half_N, prev_half_N, shif_Block, step, shift, Front, Back, temp;
//	double doun_N, current_PI_By_N, current_PI_half_By_N, step_Phase;
//	if(N :=:= 1) return;
//	double *Array:=new double[Nshl1];
//	double *Curent_Cos:=new double[Nshr1];
//	Front:=0;
//	step:=0;
//	Back:=N;
//	doun_N:=N;
//	half_N:=N;
//	j:=N;
//	step:=1;
//	prev_half_N:=N;
//	for(i:=0; i<N; i++)
//	begin
//			Array[i]:=block[i];
//	end;
//	while(j > 1)
//	begin
//		half_N:=half_N shr 1;
//		current_PI_half_By_N:=PI_half / (double)prev_half_N;
//		current_PI_By_N:=0;
//		step_Phase:=current_PI_half_By_N*2.0;
//		for(i:=0; i<half_N; i++)
//		begin
//			Curent_Cos[i]:=0.5 / cos(current_PI_By_N+current_PI_half_By_N);
//			current_PI_By_N +:= step_Phase;
//		end;
//		shif_Block:=0;
//		for(int x:=0; x<step; x++)
//		begin
//			for(i:= 0; i<half_N; i++)
//			begin
//				shift:=prev_half_N-i-1;
//				Array[Back+shif_Block+i]:=Array[Front+shif_Block+i]+Array[Front+shif_Block+shift];
//				Array[Back+shif_Block+half_N+i]:=(Array[Front+shif_Block+i]-Array[Front+shif_Block+shift])*Curent_Cos[i];
//			end;
//			shif_Block +:= prev_half_N;
//		end;
//		temp:=Front;
//		Front:=Back;
//		Back:=temp;
//		j:=j shr 1;
//		step:=step shl 1;
//		prev_half_N:=half_N;
//	end;
//	step:=step shr 2;
//	j:=j shl 1;
//	prev_half_N:=prev_half_N shl 1;
//	half_N:=prev_half_N;
//	prev_half_N:=prev_half_N shl 1;
//	while(j<N)
//	begin
//		shif_Block:=0;
//		for(int x:=0; x<step; x++)
//		begin
//			for(i:=0; i<half_N-1; i++)
//			begin
//				Array[Back+shif_Block+(i shl 1)]:=Array[Front+shif_Block+i];
//				Array[Back+shif_Block+(i shl 1)+1]:=Array[Front+shif_Block+half_N+i]+Array[Front+shif_Block+half_N+i+1];
//			end;
//			Array[Back+shif_Block+((half_N-1) shl 1)]:=Array[Front+shif_Block+half_N-1];
//			Array[Back+shif_Block+(half_N shl 1)-1]:=Array[Front+shif_Block+half_N+half_N-1];
//			shif_Block +:= prev_half_N;
//		end;
//		temp:=Front;
//		Front:=Back;
//		Back:=temp;
//		j:=j shl 1;
//		step:=step shr 1;
//		half_N:=prev_half_N;
//		prev_half_N:=prev_half_N shl 1;
//	end;
//	for(int i:=0; i<N; i++)
//	begin
//		block[i]:=Array[Front+i];
//	end;
//	block[0]:=block[0] / doun_N;
//
//	delete []Curent_Cos;
//
//	delete []Array;




//shrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshrshr>
//
double *Front, *Back, *temp; // Pointers on Front, Back, and temp massives of iteration.
double *Array:=new double[Nshl1]; // Massive with doubling size.
double *Curent_Cos:=new double[Nshr1]; // store for precomputing.

Front:=Array; // Front pointed on begin Massive.
step:=0;
Back:=Front+N; // Back pointed on center of Massive.

for(i:=0; i<N; i++)
begin
    Front[i]:=block[i]; // Load input Data
end;

while(j > 1)// Cycle of iterations Input Butterfly
begin
    half_N:=half_N shr 1;
    current_PI_half_By_N:=PI_half / (double)prev_half_N;
    current_PI_By_N:=0;
    step_Phase:=current_PI_half_By_N*2.0;
    for(i:=0; i<half_N; i++)
    begin
        //Precompute Cosine's coefficients
        Curent_Cos[i]:=0.5 / cos(current_PI_By_N+current_PI_half_By_N);
        current_PI_By_N +:= step_Phase;
    end;
    shif_Block:=0;
    for(int x:=0; x<step; x++)
    begin
        for(i:= 0; i<half_N; i++)
        begin
            shift:=shif_Block+prev_half_N-i-1;
            Back[shif_Block+i]:=Front[shif_Block+i]+Front[shift];

            Back[shif_Block+half_N+i]:=
              (Front[shif_Block+i]-Front[shift])*Curent_Cos[i];
        end;
        shif_Block +:= prev_half_N;
    end;
    temp:=Front;
    Front:=Back;
    Back:=temp;
    j:=j shr 1;
    step:=step shl 1;
    prev_half_N:=half_N;
end;

while(j<N) // Cycle of Out ButterFly
begin
    shif_Block:=0;
    for(int x:=0; x<step; x++)
    begin
        for(i:=0; i<half_N-1; i++)
        begin
            Back[shif_Block+(i shl 1)]:=Front[shif_Block+i];
            Back[shif_Block+(i shl 1)+1]:=
               Front[shif_Block+half_N+i]+Front[shif_Block+half_N+i+1];
        end;
        Back[shif_Block+((half_N-1) shl 1)]:=Front[shif_Block+half_N-1];
        Back[shif_Block+(half_N shl 1)-1]:=Front[shif_Block+(half_N shl 1)-1];
        shif_Block +:= prev_half_N;
    end;
    temp:=Front;
    Front:=Back;
    Back:=temp;
    j:=j shl 1;
    step:=step shr 1;
    half_N:=prev_half_N;
    prev_half_N:=prev_half_N shl 1;
end;
for(int i:=0; i<N; i++)
begin
    block[i]:=Front[i]; // Unload DCT coefficients
end;
block[0]:=block[0] / dN; // Compute DC.