clear
N =16;
M=N/2;
%w = .42 - .5*cos(2*pi*(0:M-1)/(M-1)) ...
%      + .08*cos(4*pi*(0:M-1)/(M-1));  %blackman window i.e. tight frame
 w=ones(1,N);  
w=w';
y = wavread('think.wav');

L=2^16;
j=zeros(1,(L/2+1000));  %adding a tone to half the signal
n=[0:(L/2+1000),j];
x=.4*sin(pi*n/8)';
x=x+y(1:2000+L+1,1);
wavplay(x,44000)



nframes=67000;


Xtwz = zeros(N,nframes);   % pre-allocate STFT output array
M = length(w);       % M = window length, N = FFT length
zp =zeros(N-M,1) ;      % zero padding (to be inserted)
R=1;                    % padding
M1=M/2;
xoff = 0;



% current offset in input signal x
for m=1:nframes
  xt = x(xoff+1:xoff+M);
  % extract frame of input data
  xtw = w .* xt;         % apply window to current frame
  xtwz = [xtw(M1+1:M); zp; xtw(1:M1)]; % windowed, z-padded
 Xtwz(:,m) = fft(xtwz); % STFT for frame m
  xoff = xoff +R;       % advance in-pointer by hop-size R
end;
colormap(hot)
%imagesc(abs(Xtwz));
subplot(411);plot(w);
subplot(411);plot(x);
subplot(412);
imagesc(abs(Xtwz(1:N,:)));
Xtwz(N/8:7*N/8,:)=0; % filter out tone
% Xtwz(N-20:N-5,:)=0; % removing other frequncys
% Xtwz(5:20,:)=0;     % removing other frequncys
subplot(413);
imagesc(abs(Xtwz(N/2:N,:)));


X=real(sum(Xtwz)')/M;
subplot(414);plot(X);
wavplay(X,44000)