Is ifft(fft(x).*fft(h)) faster or conv(x,h) ?

Dear All,

I need to find out which one is faster to obtain convolution? - Linear convolution using just 'conv'? - or Obtaining Convolution from ifft(fft)? As far as I have coded and identified using tic/toc, that Matlab performs linear convolution faster than ifft(fft). But textbooks say fft is faster.

[zhian xue]

I am also interested in this problem. Did you solved the problem?

[Alim Huseynov]

That was my task

Read the audio file ‘africa-toto.wav’ of Coursework Assignment 1 using the MATLAB command “wavread”. Call the obtained signal x[n]. Let h[n] be the 37-tap FIR filter of Coursework Assignment 1 available electronically in the file “Filter coefficients.docx”.

1) Compute the linear convolution of x[n] and h[n] using Matlab function “conv”. Measure the computation time using Matlab functions “tic” and “toc” (repeat the computation 100 times to get an accurate estimate).

(4 marks)

2) Compute the linear convolution of x[n] and h[n] using the Fast Fourier Transform and measure the computation time. Compare the results to those obtained in part 1) and comment. (5 marks)

3) Study experimentally how the results in part 1) and 2) change when the length of signal x[n] increases.

and this was conclusion

Problem 1. Audio file was read and memorized in vector as per instructions. Another vector was created for storing coefficients of 37 tap filter. Then both those vectors was convolved using Matlab Software. The operation of convolution was repeated for 100 times using 'for' loop and the average time spent for this was 0.0883 seconds. This is subject to hardware and software performance of computer used for this calculation. The code is provided in Appendix Section 1.

*******************************************************************************

Problem 2. This time Convolution operation carried out using FFT technique. The average time out of 100 operations was 1.12 seconds which was 15 times more than above operation. The code is provided in Appendix Section 2A.

The average number of operations carried out in Problem 1. was: (N1*N2)=6570432*37= 243,105,984

But in Problem 2 before applying FFT technique we carried out zero padding and then the length of both signal were 6570468. And overall number of operations were (5*N*LOG2(N))=6570468*22.6476=148,805,331=744,026,655 which was 3 times more.

The other possible reason of FFT taking long time might be the large length of signal and handling (storing) during calculations. Because on first problem computer handles two vectors with one very long and one very short. But on second case both the signals are very long. Another important point is that during usage of FFT computer handles complex numbers but on simple linear convolution involved just real numbers.

However if two other signals taken with same length N=10000 and used for obtaining convolution, it can be shown that Matlab uses less time when FFT technique is applied rather than using just normal 'conv' function. The time used by FFT was at least twice less than 'conv' function. The code in Section 2B.

But if length of one of two other arbitrary signal input kept long, say equal to 10,000 and another one incremented from lowest to 10,000, Then it can be seen that till the length of second signal becomes equal to 2000 normal convolution spends less time, and afterwards the fft becomes more efficient spending less time respectively on each increment of length of second signal. Time spent on fft technique does not change with increments of length of second signal.

The code in Section 2C.

Spikes on green (FFT Technique) are due to different algorithms applied during calculations. When the Numbers There is another type of convolution which is called Circular convolution and used normally for periodic signals. But technique can be applied for non-periodic signals if the both input signals are zero padded to some length. And this Circular convolution can be effectively calculated (replaced as the results are same) by applying FFT technique.

*************************************************************************************

Problem 3. This experiment confirms the experiment in Problem 2. Time spend on calculation increases directly proportional with increased number of elements in signal and especially for FFT. On previous one the length of signals were brought closer and closer. And while they get closer it was observed that Normal Linear convolition operation takes long time but fft technique takes less time. But in this one, the lenght of signals gets far away from each other. So the application of fft technique takes longer and longer, however normal linear convolution keeps short time value.

As seen from the plot below which characterizes average time spent with increasing number of elements. When number of elements reached 50*E5 it hit maximum, 4.5 seconds. Meanwhile Linear convolution made just in 0.1 seconds. The code is provided in Appendix A Section 3.

Sudden falls and rises within Operation Time graph of FFT technique are due to different FFT calculation algorithms chosen by Matlab. e.g. If the length of Signals can be divided into power of 2, then the algorithm of FFT called Cooley-Tukey is applied and the process accelerates considerably.

PS. The total time for overall calculation and obtaining these plot takes about 3-4 minutes. And even this duration might be longer depending on performance of computer.

***********************************************************************

[Alim Huseynov]

Section 1. (for Problem 1)

format LONGENG;          %increased precision 
[x,fs]=audioread('africa-toto.wav');
x=x';            %Column vector was transposed for making it row vector. 
yy37=[1:37];       % 37 place line vector created for further storage
h=ones(1,18);         % this line vector will accept given values
h(1)=-0.0000000000000020464110886333966;
h(2)=0.020185108207322531;
h(3)=-0.0000000000000032166461686436967;
h(4)=0.0063932717248746463;
h(5)=-0.0000000000000042068450824985665;
h(6)=-0.030584021511458139;
h(7)=-0.0000000000000037927619003410752;
h(8)=0.0092608976386982338;
h(9)=-0.0000000000000034686968012612998;
h(10)=0.055319499353126252;
h(11)=-0.0000000000000029165858917179788;
h(12)=-0.062431158774528053;
h(13)=-0.0000000000000030126051803342088;
h(14)=-0.079105057829546827;
h(15)=-0.0000000000000030966220578734095;
h(16)=0.3013776867026185;
h(17)=-0.0000000000000028385702197172922;
h(18)=0.58901882620106349;
yy37(1)=0.0053712613864042745;       %the value manually given because there is no y(0)
yy37(2:19)=h;                        %line vector copied to filter coefficients
 
for i=1:18                           % for loop used for the rest of values.
yy37(19+i)=yy37(19-i);
end
 
t = zeros(1,100);      % Vector for storing time values
 
for n = 1:100          %Looping repeats 100 times
    tic;               % Timer started
    x_c=conv(x,yy37);  %Operation of convolution
    t(n)=toc;          %During each loop relevant value of timer will be copied and timer will be reset
    
end
 
t_aver=mean(t);      %Average time spent for operation

********************************************************************************

Section 2A (Problem 2)

 format LONGENG;          %increased precision 
[x,fs]=audioread('africa-toto.wav');
x=x';            %Column vector was transposed for making it row vector.  
yy37=[1:37];       % 37 place line vector created for further storage
h=ones(1,18);         % this line vector will accept given values
h(1)=-0.0000000000000020464110886333966;
h(2)=0.020185108207322531;
h(3)=-0.0000000000000032166461686436967;
h(4)=0.0063932717248746463;
h(5)=-0.0000000000000042068450824985665;
h(6)=-0.030584021511458139;
h(7)=-0.0000000000000037927619003410752;
h(8)=0.0092608976386982338;
h(9)=-0.0000000000000034686968012612998;
h(10)=0.055319499353126252;
h(11)=-0.0000000000000029165858917179788;
h(12)=-0.062431158774528053;
h(13)=-0.0000000000000030126051803342088;
h(14)=-0.079105057829546827;
h(15)=-0.0000000000000030966220578734095;
h(16)=0.3013776867026185;
h(17)=-0.0000000000000028385702197172922;
h(18)=0.58901882620106349;
yy37(1)=0.0053712613864042745;       %the value manually given because there is no y(0)
yy37(2:19)=h;                        %line vector copied to filter coefficients
 
for i=1:18                           % for loop used for the rest of values.
yy37(19+i)=yy37(19-i);
end
  
t = zeros(1,100);      % Vector for storing time values
 
yy37Pad=[yy37 zeros(1,length(x)-1)];   %Zero Padding 
xPad=[x zeros(1,length(yy37)-1)];         %Zero Padding
 
 
for n = 1:100          %Looping repeats 100 times
    tic;               % Timer started
x_c2=ifft(fft(xPad).*fft(yy37Pad));      % Obtaining inverse transform of fft of data set.
t(n)=toc;
 end
 
t_aver=mean(t);    %Average time spent for operation

*******************************************************************************

Section 2B (Problem 2)

x = randi([1 1000],1,10000);
y = randi([1 1000],1,10000);
time1=zeros(1,20);
time2=zeros(1,20);
n = length(x) + length(y) - 1;
 
 
for i=1:10
tic
z1 = conv(x,y);
time1(i)=toc;
end
 
 
 
for i=1:10
tic
z2 = ifft(fft(x,n) .* fft(y,n)); 
time2(i)=toc;
end
 
t1_normal=mean(time1)
t2_fft=mean(time2)

*******************************************************************************

Section 2C (Problem 2)

x = randi([1 100],1,10000);
y = randi([1 100],1,10000);
 
time1=zeros(1,5);
time2=zeros(1,5);
 
time_l=zeros(1,1000);
time_2=zeros(1,1000);
 
for b=1:2000
y1=y(1:b*5);
 
 
for i=1:5
tic
z1 = conv(x,y1);
time1(i)=toc;
end
 
n = length(x) + length(y1) - 1;
for i=1:5
tic
z2 = ifft(fft(x,n) .* fft(y1,n)); 
time2(i)=toc;
end
 
time_1(b)=mean(time1);
time_2(b)=mean(time2);
end
 
figure(1)
 
plot(5*[1:2000],time_1,'r')
hold on
plot(5*[1:2000],time_2,'g')
legend('normal','fft')
xlabel('Length of second signal')
ylabel('Time spent')
set(findall(gcf,'-property','FontSize'),'FontSize',14); %size
title('Plot Showing how time spent on convolution differs with technique applied')

************************************************************************************************************************

Section 3

format LONGENG;          %increased precision 
[x,fs]=audioread('africa-toto.wav');
x=x'; % Column vector transposed
yy37=[1:37];       % 37 place line vector created for further storage
h=ones(1,18);         % this line vector will accept given values
h(1)=-0.0000000000000020464110886333966;
h(2)=0.020185108207322531;
h(3)=-0.0000000000000032166461686436967;
h(4)=0.0063932717248746463;
h(5)=-0.0000000000000042068450824985665;
h(6)=-0.030584021511458139;
h(7)=-0.0000000000000037927619003410752;
h(8)=0.0092608976386982338;
h(9)=-0.0000000000000034686968012612998;
h(10)=0.055319499353126252;
h(11)=-0.0000000000000029165858917179788;
h(12)=-0.062431158774528053;
h(13)=-0.0000000000000030126051803342088;
h(14)=-0.079105057829546827;
h(15)=-0.0000000000000030966220578734095;
h(16)=0.3013776867026185;
h(17)=-0.0000000000000028385702197172922;
h(18)=0.58901882620106349;
yy37(1)=0.0053712613864042745;       %the value manually given because there is no y(0)
yy37(2:19)=h;                        %line vector copied to filter coefficients
for i=1:18                           % for loop used for the rest of values.
yy37(19+i)=yy37(19-i);
end
% The lenght of x will be step by step increased (65 steps) and during
% each step both Linear and FFT willl be used for obtaining convolution.
% The mean time of each operation on total 3 will be stored for further
% discussion
 
M1_time=zeros(65,1);   % time variable vector for obtaining mean time during each cycle from Linear 
M2_time=zeros(65,1);   % time variable vector for obtaining mean time during each cycle from FFT
 
 
 
for k=1:65        % Loop will be repeated 65 times, each time value of k will be incremented by +1 starting from 1.
    
    x_new=x(1:100000*k); % During Each loop the lenght of x_new signal will be increased by 100000*k samples. 
        
t1 = zeros(1,3);       % Vector for storing time values
    for n = 1:3           %Looping repeats 3 times
    tic;                    % Timer started
    x_c=conv(x_new,yy37);  %Operation of convolution
    t1(n)=toc;          %During each loop relevant value of timer will be copied and timer will be reset
    end
 
t2 = zeros(1,3);      % Vector for storing time values
yy37Pad=[yy37 zeros(1,length(x_new)-1)];   %Zero Padding 
    xPad=[x_new zeros(1,length(yy37)-1)];         %Zero Padding
 
    for n = 1:3          %Looping repeats 3 times
    tic;               % Timer started
    
    x_c2=ifft(fft(xPad).*fft(yy37Pad));      % Obtaining inverse transform from fft
    t2(n)=toc;
    end
  
%Average value of time spent for operation obtained and stored in respective element of row vector. 
M1_time(k)=mean(t1); 
M2_time(k)=mean(t2);
 
end
  
figure(1)
plot(M1_time)
hold on
plot(M2_time)
xlabel('Number of elements, in 10^5')
ylabel('Time, sec')
legend ('Linear Conv','FFT')

[Rik]

Just a note: format LONGENG will not increase calculation precision, it will only affect how values are shown in the command window.

Answer 1

[Rik]

Why would a textbook say ifft(fft()) is faster? That doens't make sense. If that was the case, Mathworks would have implemented conv a bit like this:

function out=conv(x,h)
out=ifft(fft(x).*fft(h));
end

The mere fact that they didn't should tell you the real conv function is faster than the one I put here.

Answer 2

[Honglei Chen]

Time is not the best way to compare the two approaches. Rather, the best approach to describe is the computation complexity, i.e., when the size of the signal increases, how does the time of computation increases accordingly. It is in this comparison that the FFt approach shows the advantage. The theory should be available in any standard DSP text book and here is a webpage for a summary

https://ccrma.stanford.edu/~jos/ReviewFourier/FFT_Convolution_vs_Direct.html

HTH

[Rik]

I would expect the implementation in Matlab to be a black box: it is not relevant if internally there is an FFT or direct approach, as long as it gives the correct output. That is true in general for internal functions. Your comment about computational complexity still holds, but I expect the Mathworks employees that engineered the conv function to know about this.