fft_y4 = fft(y4);   % discrete Fourier transform
M4     = length(fft_y4);
w      = linspace(0, fs, M4);

figure(4); clf;
plot(w, abs(fft_y4));
grid on;
title('Spectrum of periodic signal');
xlabel('frequency (Hz)');
set(gca, 'XLim', [0 fs/2]);
