CircularConvolutionUsingDCTandDST.m

%circular convolution using Discrete sine and cosine transforms.

%Reference: V. G. Reju, S. N. Koh and I. Y. Soon, Convolution Using Discrete Sine and Cosine Transforms,
%IEEE Signal Processing Letters, VOL. 14, NO. 7, JULY 2007, pp.445-448.
%Email: [e-mail]

clear; clc;
N=20; %length of the sequences.
s=rand(N,1); %first sequence.
h=rand(N,1); %second sequence.


%Calculation of T1
T1=dC2e(s).*dC2e(h) - dS2e(s).*dS2e(h);

%calculation of T2
T2=dS2e(s).*dC2e(h) + dC2e(s).*dS2e(h);

%after calculating T2, remove zero th and Nth values from T2 because for
%S1e the range is from n=1:N-1 but T2 calculated is for n=0:N
T2=T2(2:end-1);

%Multiply the first and last element of T1 by k_{k}=2
T1(1)=T1(1)*2;  %T1(1) will have some nonzeros value, which corresponds to the DC component. So if we don's multiply , there will be a DC shift in the final convolution output. So in BSS, the scaling on this component may not cause any problem.
                % If the signals to be convolved are of zero mean, there won't be any dc and hence T1(1) will be zero.
T1(N+1)=T1(N+1)*2; %T1(N+1) is always zero, so there is no effect in multiplying but for the compatability with the std equation we can multiply.

%Calculation of dT1c1e=dC1e{T1}
% dT1c1e=dC1e*T1; 
T1C1e=DCT1e(T1);
if rem(N,2)==0  %For even N
    temp=T1C1e(1:2:end-2);
    dT1C1e=[temp; T1C1e(end); flipud(temp)];
else
    temp=T1C1e(1:2:end-1);
    dT1C1e=[temp; flipud(temp)]; 
end

%Calculation of dT2s1e=dS1e{T2}
%dT2s1e=dS1e*T2; 
T2S1e=DST1e(T2);
if rem(N,2)==0
    temp=T2S1e(2:2:end-1);
    dT2S1e=[temp; 0; -1*flipud(temp)];
else
    temp=T2S1e(2:2:end);
    dT2S1e=[temp; -1*flipud(temp)];
end


y=1/(8*N)*(dT1C1e(2:end)+[dT2S1e;0]);

%Convolution by DFT method
ydft= real(ifft(fft(s).*fft(h)));

figure(1);
%DTT method
plot(y); hold on;
%DFT method
plot(ydft,'r*');hold off;
legend('Proposed DTT method','DFT method');
title('Convolution output using DFT method and proposed DTT method, to verify the proposed method.');