# Plot the sum frequency generation spectrum using convolution MATLAB

I am attempting to calculate the spectrum of a pulse that has undergone sum-frequency generation (in this case it is a gaussian, so it is correct to also say frequency doubling/Second harmonic generation). The SHG signal in the frequency domain is given as,

$$E_{SHG}(2\omega) = E_1(\omega)*E_2(\omega)$$

Therefore a signals SHG spectrum is just an autoconvolution of the original spectrum.

However, I am unfamiliar with the practical use of discrete convolution and do not know how to transform the new x-axis in to a suitable vector for plotting?

clear all; close all;
dt = 0.01;
x = 200:dt:1000;    %Frequency axis (THz)

%Generate Stokes Profile
width_stokes = 20;

% center frequency
f = 500;
Es = exp(-(x-f).^2/width_stokes^2);
Es=Es./max(Es);
plot(x,Es);
title("Stokes spectrum");

SHG = conv(Es,Es,'same');
SHG = SHG./max(SHG);

figure
% New x-axis for SHG plot
x1 = (1:length(SHG));
plot(x1,SHG)

xlabel('frequency (A.U.)')
$$$$


You have a gaussian centered at 500 THz. We would expect the convolution to have a single gaussian centered at 1000 THz.

A linear convolution of two sequences of N points each will have a length of 2*N-1 samples. You have the added complication that your frequency vectors don't start at 0Hz. One way to fix this would be to have them simply start at zero

dt = 0.01;
x = 0:dt:1000;    %Frequency axis (THz) starting at 0
N = length(x);
freqAxisConvolution = dt*(0:2*N-2);


Alternatively, you can just calculate the offset. If both vectors where unit impulses (starting at 200 THz) the convolution would be a unit impulse (starting at 400 THz), so the offset is simply the sum of the individual offsets. In other words if each vectors spans from 200 THz to 2000 THz the convolution will span from 400 THz to 4000 THz

Here is the full thing

%%
close all
dt = 0.01;
x = 200:dt:1000;    %Frequency axis (THz)

%Generate Stokes Profile
width_stokes = 20;

% center frequency
f = 500;
Es = exp(-(x-f).^2/width_stokes^2);
Es=Es./max(Es);
plot(x,Es);
title("Stokes spectrum");

% discrete convolution produces 2*N-1 output samples
SHG = conv(Es,Es,'full');
SHG = SHG./max(SHG);

figure
%  X axis: spans the sum of the original axes
N = length(Es);

freqAxis = 2*x(1)+(0:2*N-2)*dt;

plot(freqAxis,SHG)
grid on
xlabel('frequency in THz');
`
• Thank you so much for the help. The explanation makes sense and works! The main problem for me was not starting from zero in the frequency array, thanks again.
– Ryan
Commented Jan 12, 2022 at 18:23