M = 200;
L=M+1;
n=[0:M]-M/2;
D = 8;
fc=1/(2*D); h=(2*fc)*sinc(2*fc*n);
NFFT = 2^12; f=[0:NFFT-1]/NFFT;

figure(1);
subplot(211); plot(f,abs(fft(h,NFFT)));
h = h.*hamming(L).';
subplot(212); plot(f,abs(fft(h,NFFT)));

B = diag(h)*exp(1j*2*pi*[0:M].'*[0:D-1]/D);
figure(2);
subplot(211); plot(f,abs(fft(B,NFFT)));
subplot(212); plot(f,20*log10(abs(fft(B,NFFT))));
subplot(211);
hold on;
ax = axis;
for i=1:D
    m = mod(-(i-1),D);
    text(m/D,0.5,int2str(i-1),'FontSize',16);
end
hold off;

P = ceil(L/D);
L = P*D;
g = zeros(1,L);
g(1:length(h)) = h;
G = reshape(g,[D,P]);

if 0
    [x,Fs]=audioread('fireflyintro.wav');
else
    Fs = 40e3;
    t = [0:2*Fs];
    fi = t/(2*Fs);
    phi = 2*pi*cumsum(fi);
    x = exp(1j*phi);
end
T = D*floor(length(x)/D);
X = reshape(x(1:T),[D,T/D]);
X = flipud(X);

Y = zeros(size(X));
for i=1:D
    Y(i,:) = conv(X(i,:),G(i,:),'same');
end

Y = fft(Y,D);

figure(4);
spectrogram(x,hamming(512*4),485*4,1024*4,Fs,'yaxis');
hold on;
ax = axis;
for i=1:D
    plot(ax(1:2),((i-1)/D + 1/(2*D))*(Fs/1000)*[1 1],'r');
    m = mod(-(i-1),D);
    text(1,(m/D)*(Fs/1000),int2str(i-1),'FontSize',16);
end
hold off;

figure(3);
for k = 0:D-1 % Zero based
    spectrogram(Y(k+1,:),hamming(512),485,1024,Fs/D,'yaxis');
    title(int2str(k),'FontSize',16);
    pause;
end