Constant Multiple in trapz() result?

I'm trying to calculate an integral numerically using trapz function. However, my result is off from the expected result by a constant multiple. I was wondering if there's something about the trapz function that I dont understand. More details:

I'm trying to plot the function V and also calculate its Fourier Coefficients and plot them given the following formulae. Note that I use the cosine and sine series for Fourier series instead of the exponential series in the fomula below. Also, I agnore the multiples outside the squareroot in both the original and reconstructed functions.

I call the following code to calculate the integral and assigne the values to a_n and b_n for, say, 1000 values of n. I've added some comments to make the code clearer without the detail of the problem context. The result that I get in VSWR is not a reconstruction of the above signal as intended. In addition, my Fourier coefficients as calculated are not getting smaller with increasing 'n' as expected. Can you say if the use of trapz function is the problem and if the implementation results in some offset due to the loop variable n being a variable in the integral? Something like the following point in the trapz documentation?

copied from: https://www.mathworks.com/help/matlab/ref/trapz.html

• If X is a scalar spacing, then trapz(X,Y) is equivalent to X*trapz(Y).

You'll notice that the following code does the following:

1. Generate a swept frequency cosine waveform from start_freq to end_freq with the default number of sampels.
2. Resample this waveform to create an upsampled waveform that is the same size as the f_range vector. This ensures that the sampled waveform.
3. The waveform is used to calculate Fourier coefficients a_n, b_n.
4. Reconstructed waveform is synthesized using the Fourier coefficients.

Note that the resultant waveform has the correct frequency as the original swept frequency waveform (if the number of Fourier coefficients that is calculated is large enough) but it is biased above the original waveform by a constant factor. I have tried to use various values of num_FourierCoefficients and that changes the bias / offset magnitude.

Also note that the Fourier coefficients should be monotonically decreasing in magnitude but that does not happen in our case here.

%Calculation of Fourier Series Coefficients for the representative function
%of square of VSWR. Refer to the documentation file for mathematics


%Calculation is done based on a formula developed in reference file and
%calculates first 'num_FourierCoeff' fourier coefficients of the squared VSWR
%function which is then reconstructured the non-zero coefficients among
%these. 

%Constants
c = 2.998e8; % (m/s) - speed of light
eta = 120*pi;


%Used values of frequency and path distance for this calculation. These
%values may change for specific cases. However, we are trying to understand
% how the frequency component from the n-th harmonic --DC, first harmonic,
% second harmonic and so on vary with 'n' and how large a value of n is
% enough to capture VSWR with required precision.
h_T = 40; % (m) - Transmitter drone altitude
h_R = 40; % (m) - Receiver drone altitude
d_S = 50; % (m) - Horizontal distance between the two drones
R_d = sqrt(d_S^2 + (h_T - h_R)*(h_R - h_T)); %(m) - the line of sight distance between the transmitter and receiver (direct path)
theta_i = atan(d_S/(h_T + h_R));              % (rad.) - the angle of incidence also equal to angle of reflection
R_TS = h_T/cos(theta_i); % (m) - the distance from the transmitter to point of reflection on soil
R_SR = h_R/cos(theta_i); % (m) - the distance from point of reflection on soil to receiver

delta_R = R_TS + R_SR - R_d; %(m) - the path difference between direct and reflected path
% In the following line, Set to 1 since it doesn't effect the frequency
% component of harmonics, only their amplitude.
P_T = 1e-3; % (W) - Transmitted power = 1 mW

%Power delivered to the antenna for radiation (by Balanis, page 82, eq.
%2-77) is: P_T = 1/2*V^2/Z_0 where we use the impedance of free space as
%Z_o for the following case. Note that I'm assuming 100% radiation
%efficiency such that all the power delivered to the antenna is radiated
%into free space (air) and is used for our radar transmission

V_D = sqrt(2*eta*P_T);  %(-) - magnitude of direct voltage wave
delta_f = c/delta_R;    % - period of function in frequency
%f = 30e6;               % frequency for calculation
%%%Note: We will have to incorporate the frequency dependence of the
%%%reflectivity from our previous model. That should attenuate the
%%%amplitude of the VSWR waveform but does not effect sampling frequency.
r_alpha = 0.5;          % reflectivity of the soil surface. Effects amplitude of some harmonics but not frequency.
Lambda = (R_d/(R_TS + R_SR))*sqrt(r_alpha); % reflection coefficient (attenuation) of the reflected wave 

num_FourierCoeff = 1000; % Number of Fourier Coefficients to be calculated
start_freq = 10e6;
end_freq = 24e6;
n_points = 10*end_freq/start_freq;
%f_range = linspace(start_freq, end_freq, n_points); %frequency range to calculate the VSWR
num_periods = floor((end_freq-start_freq)/delta_f); %number of periods of VSWR function in the frequency range
end_freq =start_freq+(num_periods*delta_f);
f_range = linspace(start_freq, end_freq, n_points); %Rounding down to an integer number of periods

%Fourier Coefficients a's and b's
a_n = zeros(1,num_FourierCoeff);
b_n = zeros(1,num_FourierCoeff);

%The VSWR reconstructed function using Fourier harmonics
VSWR_reconstructed = zeros(1, size(f_range,2)); 
harmonics_a_n = zeros(1, size(f_range,2)); 
harmonics_b_n = zeros(1, size(f_range,2));
%The mathematical function describing VSWR
VSWR_sweptFreq_input = frest.Chirp('Amplitude',1, 'FreqRange',[start_freq end_freq], 'FreqUnits', 'Hz', 'InitialPhase', 0);
ts1 = generateTimeseries(VSWR_sweptFreq_input);
output = resample(ts1.data,n_points,size(ts1.data,1));
func_VSWR =sqrt(1+((2*Lambda/(1+Lambda^2))*output

%Plots to make sure the sweptfrequency waveform is sampled correctly
figure(1);
plot(ts1.data)
figure(2);
plot(f_range,func_VSWR);


%VSWR_sweptFreq_input2=chirp(f_range,  start_freq/delta_f ,end_freq, end_freq/delta_f , 'linear');

%Fourier coefficients from a_1, b_1 onwards to higher harmonics
for n = 2:1:num_FourierCoeff
    a_n(n) = (1/(delta_f))*(trapz(f_range, func_VSWR.*cos(2*pi*(1/delta_f)*f_range*(n-1))))/(num_periods);
    b_n(n) = (1/(delta_f))*(trapz(f_range, func_VSWR.*sin(2*pi*(1/delta_f)*f_range*(n-1))))/(num_periods);
end
%Fourier coefficients a_0, b_0
a_n(1)=(1/(delta_f*2))*(trapz(f_range,func_VSWR))/num_periods; 
b_n(1)=0;
%set really small values to zero 
a_n(abs(a_n)<1e-4) = 0;
b_n(abs(b_n)<1e-4) = 0;

for m = 1:1:num_FourierCoeff
    harmonics_a_n = a_n(m) * cos(2*pi*(delta_R/c)*f_range*(m-1));
    harmonics_b_n = b_n(m) * sin(2*pi*(delta_R/c)*f_range*(m-1));
    VSWR_reconstructed = VSWR_reconstructed +  harmonics_a_n + harmonics_b_n;
end

%Save Fourier coefficients in a file for later analysis
fileID = fopen('C:\Farah Akhtar\Bistatic Radar Active Soil Moisture Sensing -BRASS\Fourier Series Calcuation (Matlab)\FourierCoefficents_f10MHz_v3.txt','w');
for n = 1:1:size(a_n)
    fprintf(fileID, 'n = %1 a_n = %0.3f, b_n = %0.3f\n', n, a_n(n), b_n(n));
end
fclose(fileID);

% Rendering the Original VSWR mathematical function plotted along with the
% reconstructed VSWR on the same axis
figure(3);
plot(f_range,func_VSWR,'r',f_range,VSWR_reconstructed,'b--');
title('Original and Reconstructed VSWR functions')
xlabel(['Frequency' newline 'Range:' num2str(start_freq/1e6) ' to ' num2str(end_freq/1e6) ' (MHz)' ])
ylabel('Voltage (V)')
legend('Original','Reconstructed')

size(a_n)
size(b_n)

% plot only non-zero Fourier Coefficients
if(size(nonzeros(a_n),1) > size(nonzeros(b_n),1))
    a_n = nonzeros(a_n);
    b_n = b_n(1:size(a_n))';
else
    b_n = nonzeros(b_n);
    a_n = a_n(1:size(b_n))';
end

figure(4);
subplot(2,1,1);
stem(1:size(a_n,1), a_n);
title('Fourier Cofficients a''s')
xlabel('n')
ylabel('Magnitude')
subplot(2,1,2);
stem(1:size(b_n,1), b_n);
title('Fourier Cofficients b''s')
xlabel('n')
ylabel('Magnitude')

[Ameer Hamza]

Several variables are not defined in this code: delta_f, delta_R, Lambda. Without running this code with correct parameters, it is difficult to suggest anything.

[Farah Akhtar]

Thanks for your response, Ameer Hamza. I was trying to avoid complicating the question because of how it is a part of a larger problem. I have updated the code again to illustrate what I was trying to do.