h(1:128)=cos((0:127)*2*pi*6/256);
h(31:71)=h(31:71)+0.5*cos((-10:30)*2*pi*52/256);
h(129:256)=cos((0:127)*2*pi*25/256);
N = 256;
H = [fft(h) fft(h) fft(h)];
al = [0:(N/2)-1, -(N/2):-1];
for n=1:N
W=exp(-2*pi^2*al.^2./((n-(N/2)-1).^2));
S(n,:) = ifft(H(n+(N/2):n+(N/2)+N-1) .* W);
end