pitch_shift.m

load "brain1_denoised.mat"

% test signals
%
%result = zeros(1, 512);
%result = sin(2*pi*(1:2048)*20/512) + cos(2*pi*(1:2048)*6/512);

result = (result-mean(result))/max(result);

printf("Denoised brain data loaded\n");
max(max(result))

srate = 512;
sound_srate = 44100;

eeg = result;
sl = min(1048576, length(eeg)); % 512 gives you a second
%sl = min(2048, length(eeg)); % 512 gives you a second

ws = ceil(1.0*srate); % 0.5 seconds
sound_max_ws = ceil(ws*sound_srate/srate)

scaling = 128;

hanning = 0.5*(1-cos(2*pi*(0:(ws-1))/(ws-1)));
sound_hanning = 0.5*(1-cos(2*pi*(0:(sound_max_ws-1))/(sound_max_ws-1)));

n_window_phases = 24
wd = round(ws/n_window_phases); % window difference
hanning_sum = 0.5*n_window_phases;

result = zeros(1, ceil(sl + 2*ws+2)*sound_srate/srate);
size(result);

eeg_start = 1-wd;

window_starts = [1 + round(ws*(-1:(n_window_phases-2))/ n_window_phases)];
wc = 0;

printf('Precomputing sound waves\n');

for i=1:(ws/2)
  cplexp_fundamental_lookup = exp(j*2*pi*(0:(sound_max_ws-1))./sound_max_ws);
end

size(cplexp_fundamental_lookup)
sound_idx_vector = 0:(sound_max_ws-1);
size(sound_idx_vector)
%hannings = zeros(50, length(result));
while (min(window_starts) < sl)

  for p=1:n_window_phases
    % window phases
    window_eeg_start_index = window_starts(p) - eeg_start + 1;
    window_eeg_end_index   = window_starts(p) - eeg_start + ws;
    window_time_start = round((-eeg_start + window_starts(p))*sound_srate/srate + 1);

    if (window_eeg_end_index > length(eeg))
      % Go off the end
      signal = eeg(window_eeg_start_index:end);
      signal = postpad(signal, ws);
    elseif (window_starts(p) < 1)
      signal = eeg(1:window_eeg_end_index);
      signal = prepad(signal, ws);
    else
      signal = eeg(window_eeg_start_index:window_eeg_end_index);
    endif

    if (prod(size(signal)) == 0)
      signal = zeros(1, ws);
    end
    signal = signal.*hanning;

    % Now perform the pitch shifting.
    signal_fft = fft(signal);
    %signal_fft

    printf('fft evaluated\n');
    % resample at higher pitch
    window_result = zeros(1, sound_max_ws);

    for si = 0:(length(signal_fft)/2-1)
      %
      cplexps = 2*real(cplexp_fundamental_lookup(1+mod(si*sound_idx_vector*scaling+window_time_start, sound_max_ws)));
      window_result = window_result + cplexps.*signal_fft(si+1);
    endfor

    printf('Resampled\n');

    % size(window_result)

    window_result = real(window_result).*sound_hanning;
    n = 1:sound_max_ws;

    idx = round((-eeg_start + window_starts(p))*sound_srate/srate + n);
%    printf('');
    result(idx) = result(idx) + window_result./hanning_sum;
%    hannings(n_window_phases*wc+p, idx) = window_result;

%    for n=1:sound_max_ws
%      idx;
%%      hannings(3*wc+p, idx) = sound_hanning(n);
%      result(idx) = result(idx) +
%      real(window_result(n))/
%%      hannings(3*wc+p, idx) =     real(window_result(n))/(sound_hanning(1+mod(n-1+round(window_starts(1)*sound_srate/srate), sound_max_ws)) + sound_hanning(1+mod(n-1+round(window_starts(2)*sound_srate/srate), sound_max_ws)) + sound_hanning(1+mod(n-1+round(window_starts(3)*sound_srate/srate), sound_max_ws)));
%    endfor

    printf('Rebuilt\n');
  endfor

  window_starts = window_starts .+ ws
  wc=wc+1;
  printf('window %d\n', wc);
  max_window_starts = max(window_starts);
  max_window_starts
  sl;
end

result;

plot(result);