- The convolution by for-loop needs to flip the kernel (Morlet)
- The MATLAB convolution 'same' dimension applied for the first array, so you need to swap arguments

# Why result of conv and point wise multiplication has different result?

8 views (last 30 days)

Show older comments

I am comparing conf function and sum over point wise multiplication of two signal and then plot it. I end up with figure that have different axis length. Literature says that total point convolution is Signal+Kernel-1, but when i check dimension of con_auto is 1x4001. Does it should be

length of cmw (4001)+length of chirpTS (6001)-1?

Does anyone can explain to me?

here is modification code that I use

clc,close; clear

srate=1000; f=[2 8];

t=0:1/srate:6; n=length(t);

% create chirp signal

chirpTS =sin(2*pi.*linspace(f(1),f(2)*mean(f)/f(2),n).*t);

% create complex Morlet wavelet

% wtime is wavelet time, cmw is complex wavelet waveform

wtime=-2:1/srate:2;

w=2*(4/(2*pi*5))^2;

cmw=exp(1i*2*pi*5.*wtime).* exp((-wtime.^2)/w);

% half of the length of the wavelet

halfwavL = floor(length(wtime)/2);

% zero-pad chir --> zero padding dilakukan agar hasil konvolusi mempunyai panjang yang sama

% dengan signal

chirpTSpad = [zeros(1,halfwavL) chirpTS zeros(1,halfwavL)];

% CONVOLUTION

convres =zeros(size(chirpTSpad));

for i=halfwavL+1:length(chirpTS)+halfwavL-1

%at each time point,compute dot product

convres(i)=sum(chirpTSpad(i-halfwavL:i+halfwavL).*cmw);

end

% cut off edges

convres =convres(halfwavL:end-halfwavL-1);

subplot(511),plot(t,chirpTS),title ('chirp signal')

subplot(512), plot(wtime,abs(cmw)), title ('morlet wavelet (abs)')

subplot(513),plot(t,abs(convres)), title ('convolution result')

%CALCULATE CONVOLUTION USING CONV Matlab Function

% 'same' - returns the central part of the convolution that is the same size as A.

con_auto=conv(cmw,chirpTS,'same');

subplot(514), plot(abs(con_auto)),

title ('convolution result using conv function')

% CONVOLUTION VIA FREQ DOMAIN

Lconv=length(t)+length(wtime)-1;

convres2=ifft( fft(chirpTS,Lconv).*fft(cmw,Lconv) );

convres2=convres2(halfwavL:end-halfwavL-1);

%line 29 also can be written as

subplot(515), plot(t,abs(convres2))

title ('convolution result via frq domain')

##### 0 Comments

### Accepted Answer

Bruno Luong
on 23 Aug 2023

Edited: Bruno Luong
on 24 Aug 2023

You have at least 2 errors in your codes

I fix both and here is the result

clc,close; clear

srate=1000; f=[2 8];

t=0:1/srate:6; n=length(t);

% create chirp signal

chirpTS =sin(2*pi.*linspace(f(1),f(2)*mean(f)/f(2),n).*t);

% create complex Morlet wavelet

% wtime is wavelet time, cmw is complex wavelet waveform

wtime=-2:1/srate:2;

w=2*(4/(2*pi*5))^2;

cmw=exp(1i*2*pi*5.*wtime).* exp((-wtime.^2)/w);

% half of the length of the wavelet

halfwavL = floor(length(wtime)/2);

% zero-pad chir --> zero padding dilakukan agar hasil konvolusi mempunyai panjang yang sama

% dengan signal

chirpTSpad = [zeros(1,halfwavL) chirpTS zeros(1,halfwavL)];

% CONVOLUTION

convres =zeros(size(chirpTSpad));

for i=halfwavL+1:length(chirpTS)+halfwavL-1

%at each time point,compute dot product

convres(i)=sum(chirpTSpad(i-halfwavL:i+halfwavL).*flip(cmw)); % BUG here, flip cmw

end

% cut off edges

convres =convres(halfwavL:end-halfwavL-1);

subplot(511), plot(t,chirpTS),title ('chirp signal')

subplot(512), plot(wtime,abs(cmw)), title ('morlet wavelet (abs)')

subplot(513), plot(t,abs(convres)), title ('convolution result')

%CALCULATE CONVOLUTION USING CONV Matlab Function

% 'same' - returns the central part of the convolution that is the same size as A.

con_auto=conv(chirpTS,cmw,'same'); % % BUG here, swap signals

subplot(514), plot(t,abs(con_auto)),

title ('convolution result using conv function')

% CONVOLUTION VIA FREQ DOMAIN

Lconv=length(t)+length(wtime)-1;

convres2=ifft( fft(chirpTS,Lconv).*fft(cmw,Lconv) );

convres2=convres2(halfwavL:end-halfwavL-1);

%line 29 also can be written as

subplot(515), plot(t,abs(convres2))

title ('convolution result via frq domain')

##### 0 Comments

### More Answers (1)

Matt J
on 23 Aug 2023

In one case, you are using conv(___'same'). The 'same' flag will throw away some of the resulting points.

In the second case, you are doing fft-based convolution, which implements circular convolution, rather than linear convolution. The relation Signal+Kernel-1 only applies to linear convolution.

##### 6 Comments

Bruno Luong
on 23 Aug 2023

Edited: Bruno Luong
on 23 Aug 2023

I'm lost. You cite a document that is not MATLAB related, or do I miss something?

Just for a fun of OT discussion?

The question is about the difference of the 4th plot; which is obviously not due to small arithmetic errors among different implementations.

Walter Roberson
on 23 Aug 2023

There is nothing anywhere in the documentation that restricts MATLAB to produce bitwise identical results to what naive MATLAB-level code would produce for the operation -- and there is a long history of cases where actual implementations produce different results than naive implementations. Actual implementations are permitted to use whatever hardware instructions are appropriate to the situation -- and considering the nature of the Execution Engine in internally producing machine code for execution, the Execution Engine is permitted to produce different hardware instruction sequences for different CPU models.

At the moment, the only restrictions that I can recall as being documented (sometimes documented only in bug reports) are:

- MATLAB restricts itself to not chaining Extended Precision float operations (such as Intel 80 bit registers) -- that it will store back the results of each operation and load again if appropriate or otherwise inject instructions so that each calculation is as-if on 64 bit registers. This restriction might possibly not apply to external libraries it calls such as LAPACK. I do not recall at the moment whether the implication is that MATLAB refrains from using FMA (Fused Multiply and Add) instructions; I think Simulink has a configuration option about whether FMA can be used by Simulink
- When additions are partitioned across multiple cores, MATLAB restricts itself to adding the partial results in a stable order, rather than as the results are available

But, for example, the Windows recommended hardware is "Any Intel or AMD x86–64 processor with four logical cores and AVX2 instruction set support" and at least as of R2023a, notes that at some point in the future AVX2 will be required -- implying that AVX2 is already being used by MATLAB when available, and thus that generally MATLAB gives itself permission to use efficient instructions if the CPU supports them.

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!