Confusion in FFT of Resampled Data I have the following code in Matlab:

fs = 1e6;

y = randn(1,100);
y = upsample(y,8);
windowSize = 8;
b = (1/windowSize)*ones(1,windowSize);
a = 1;
y = filter(b,a,y);

fft1 = fft(y);
fft1 = fft1(1:length(y)/2+1);
psd = (1/(fs*length(y)))*abs(fft1).^2;
psd(2:end-1) = 2*psd(2:end-1);
f = 0:fs/(length(y)):fs/2;
figure(1);
plot(f,10*log10(psd))

y = resample(y,7,4);
fft1 = fft(y);
fft1 = fft1(1:length(y)/2+1);
psd = (1/(fs*length(y)))*abs(fft1).^2;
psd(2:end-1) = 2*psd(2:end-1);
f = 0:7*fs/(4*length(y)):7*fs/8;
figure(2);
plot(f,10*log10(psd))


In this code I have a random noise signal and I am upsampling it and passing it through a zero order hold filter then plotting its Power spectral Density (PSD) in figure 1.
After resampling the signal to a new rate when I check its PSD in figure 2 I get some different result.
Could someone explain why is the spectrum of the resampled data different from the original one and what can I do to get the same PSD as the original signal?
As you know you have to use an ideal low-pass filter, and zero order hold is an approximation which needs more compensation. Ideal LPF has a flat response with sudden cut-off and zero stop band but your moving average window has a sinc frequency response. So the deep valleys belong to the zeros of sinc function and the remaining which has high attenuation passed from the stopband of your filter (which is not zero because it's not an ideal LPF).