function Y = freqshift(X,SR,F,phi)
% Y = freqshift(X,SR,F,phi)
%    Apply a linear frequency shift of F Hz (positive = upwards) to X.
% 2011-04-22 Dan Ellis dpwe@ee.columbia.edu

if nargin < 4;  phi = 0;  end

Wsec = 1.0;
W = 2^round(log(Wsec*SR)/log(2));

% Force W even
W = W + rem(W,2);
% Hop between windows - 50% overlap
H = W/2;

% Build 50% overlapped, windowed windows
%Xw = diag(hanning(W)')*frame(X,W,H);  % very very slow
% prepad by a half window to avoid edge effects
Xf = frame([zeros(H,1);X;zeros(H,1)],W,H);
nc = size(Xf,2);
Xw = repmat(hanning(W),1,nc).*Xf;

% Take FFT
SX = fft(Xw);
% Make analytic
SX = repmat([1:W]' < (W/2),1,nc).*SX;
% For shifting down, also remove low frequencies
if F < 0
  SX = repmat([0:W-1]' > (-F)/SR*W,1,nc).*SX;
end
% IFFT
Ywc = ifft(SX);

% Overlap-add
Yc = ola(Ywc,H);

% Apply frequency shift, take real part, compensate for lost -ve energy
Y = 2*real(Yc.*exp(j*F/SR*2*pi*[0:length(Yc)-1]'+phi));

% make sure it's long enough & trim off added padding
Y(length(X)+1+H) = 0;
Y = Y(H+[1:length(X)]);
