Segment a Gaussian into sections, then create Gaussians of conserved area in each section

My basic goal is as follows: I have some simulated data, for example a broad Guassian function with white noise (testp). I would like to take this function, segment into [1:400] section. In each section, I want to find the area, and then create a new, narrower Gaussian function in each interval.

Two problems I am having:

1.) i am using trapz to find the area, and then simulating narrower gaussians with the area found by trapz. However, it doesn't appear to be the same (new function space has about half the area as the original data)

2.) I want to keep the noise of the original data. Multiplying by a guassian, everywhere than the new function is zero the noise goes way (duh). Any suggestion/ideas for what I can add to my noise such that the orignal noise is preserved? I attached a picture of what I have so far (area not quite right, noise not preserved) and my code so far. Red is new desired function, blue is original noisy data. Also tips for cleaning my code up are welcome.

tsep_1d = [0.001:.005:30];
collfreq = 1./(tsep_1d(2) - tsep_1d(1));
PM1 = 2; dp_PM1 = PM1.*collfreq;


index = 1;

peakmz_array = zeros(dp_PM1,length(tsep_1d)./dp_PM1);
peakmz_array_mod = zeros(dp_PM1,length(tsep_1d)./dp_PM1);
for i = 1:length(tsep_1d)./dp_PM1; %index tsep for a given peak
peakmz_section = testp(index:index+dp_PM1-1); %index a given peak section
peakmz_array(:,i) = peakmz_section; %add each section to peakmz_array
%area_peak = trapz(peakmz_section); %find area in a given section
area_peak =trapz(tsep_1d(index:index+dp_PM1-1),testp(index:index+dp_PM1-1));
tpeak = [1:1:400]; %time interval per 2d modulation
mod_signal_mz = peak_gauss(1,tpeak,50,50); %modulation peak with random wb and tr
peakmz_section_mod = peakmz_section.*mod_signal_mz; %multiply each section my modulation signal
noise_peak_mz = std(testp(1:20)); %find noise per mz

peakmz_array_mod(:,i) = peakmz_section_mod;  % add each modulated section to peakmz_array_mod
index = index+400-1;
end

peakmz_reconstructed = reshape(peakmz_array,[1,size(peakmz_array,1).*size(peakmz_array,2)]) ;
peakmz_mod_reconstructed = reshape(peakmz_array_mod,[1,size(peakmz_array_mod,1).*size(peakmz_array_mod,2)]) ;



figure;
plot(tsep_1d,peakmz_reconstructed,'linewidth',2);
hold on;
plot(tsep_1d,peakmz_mod_reconstructed,'linewidth',2);

function [peaknd] = peak_gauss(A,x,wb,tr)
%creates a guassian peak o
%A area (constant)
%x is a vector of data points 
%wb = peak width in data points
%tr = retention time in data points

%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here
sigma = wb ./4;
peaknd = A*(1/sqrt(2*pi*sigma.^2) )*exp(-( (x-tr).^2  ./  (2*sigma.^2)    ) );

end