function [mag,f]=mag_cal(y,u,fs)
%
% This function computes the magnitude of the frequency response between u (input)
% and y (output) using the spectral analysis.
%
NFFT=1024;
[Py,f]=psd(y,NFFT,fs,NFFT,NFFT/2);
Pu=psd(u,NFFT,fs,NFFT,NFFT/2);
Pyu=Py./Pu;
mag=10*log10(Pyu);
plot(f,mag)