why does my integrator perform noise shaping?

Question

I am designing (in matlab) a simple system with a DC signal at 0 Hz and some added noise. I am performing a differentiation of the signal (by simply subtracting the last values from each other), next i am applying rolling average algorithm, to filter the high frequencies, and then i simply integrate over a number of samples. Using the pwelch() function on the output i noticed that there is some noise shaping going on, and i don't really now where it is coming from...

I am aware of delta-sigma modulators, but this is not a thing i am doing here: for my integrator i am using a very simple algorithm that adds $n$ last samples together. I would expect a simple integrator behaviour here (with frequency response falling 20db/dec).

Why does the noise shaping take place? https://i.sstatic.net/dTfLl.png

I am posting the code for the signal generation, differentiator and integrator below.

clear
home
close all

%generate the signal
    fs = 1;
    T = 1/fs;
    t = 0:T:(2^14);
    pos = 3 * t; 
    noise = 2 * randn(1, length(t));
    pos = pos+noise;


    %differentiate the signal
    diff_time = 1;
    diff_out = [];
    len = length(pos);
    for j = 1+diff_time : len
        diff_out(j-diff_time) = pos(j)-pos(j-diff_time);
    end

    %filter the hf with MA filter with 16 taps
    MA_avg = zeros(1,16);
    for i = 1 : length(diff_out)-16
       avg_val = MA_avg(end);
       old_val = diff_out(i);
       current_val = diff_out(i+16);
       MA_avg(end+1) = current_val-old_val+avg_val;
    end

    %integrate
    pos_intgr = integrate(MA_avg, 16);

    %get the spectral power
    NFFT = length(pos_intgr);
    [P, F] = pwelch(pos_intgr,ones(NFFT,1),0,NFFT,fs,'power');
    PdBW = 10*log10(P);
    plot(F,PdBW)
    title("pwelch")
    xlabel('Frequency')
    ylabel('Power spectrum (dbW)')

    function [acu_down] = integrate(signal, taps)
       len = length(signal);
       acu_down = [0];
       i = 1;
       k = 1;
       while i <= len-taps
          sum = 0; %reset the sum
          for j = 0 : taps-1
               sum = sum + signal(i+j); 
          end
          i = i+taps;
          acu_down(k) = sum;
          k = k+1;
      end
      %cut out the LSB
      acu_down = acu_down ./ taps;
   end

Answer 1

The dominant effect you see is from the differentiation which is a high pass function.

Observe the spectrum after differentiation alone:

As would be clearer on a log log plot the signal is going up at rate f versus frequency. (consistent with a differentiation).

Next observe the spectrum at the moving average output:

This is consistent with expected: The moving average approximates a Sinc filter response, and the nulls of a Sinc function go down at rate f versus frequency. With the differentiator going up at rate f and the moving average going down at rate f we see the flat response as observed with the nulls at $1/T$ where $T$ is the duration of the moving average.

Your final operation that you call an integrator is actually another average combined with a down-sample (decimation), and with that would result in the low frequency portion of your signal after another 16 point average; so another Sinc filter response of similar form and then the decimated (close in) spectrum out to half of the first main lobe: