function rayleigh_t=Rayleigh_filtered(fmax,fs,Ns)
N=8;
while(N)
    if(N <2 * fmax*Ns/fs)
        N=2*N;
    else
        break;
    end
end

N_inv=ceil(N*fs/(2*fmax))

f_interval=2*fmax/N;
T_inv=1/fs;

I_input_time=randn(1,N);
Q_input_time=randn(1,N);

I_input_freq=fft(I_input_time);
Q_input_freq=fft(Q_input_time);

SEZ(1)=1.5/(pi*fmax);
for j=2:N/2
    f(j)=(j-1)*f_interval;
    SEZ(j)=1.5/(pi*fmax*sqrt(1-(f(j)/fmax)^2));
    SEZ(N-j+2)=SEZ(j);
end

k=3;
p=polyfit(f(N/2-k:N/2),SEZ(N/2-k:N/2),k);
SEZ(N/2+1)=polyval(p,f(N/2)+f_interval);

I_output_freq=I_input_freq.*sqrt(SEZ);
Q_output_freq=Q_input_freq.*sqrt(SEZ);

I_temp=[I_output_freq(1:N/2) zeros(1,N_inv-N) I_output_freq(N/2+1:N)];
I_output_time=ifft(I_temp);
Q_temp=[Q_output_freq(1:N/2) zeros(1,N_inv-N) Q_output_freq(N/2+1:N)];
Q_output_time=ifft(Q_temp);

channel=I_output_time+i*Q_output_time;
channel=channel(N:Ns);

rms=sqrt(norm(channel,'fro'));
channel_unit=channel/rms;
rayleigh_t=diag(channel_unit'*channel_unit)';