%%% Real-Time Speech Pitch Shifting on an FPGA
%%% Estephan, Sawyer, Wanninger
%%% M-File by Scott Sawyer (2006)
%%% Villanova University
%%% Department of Electrical and Computer Engineering


clear; clc;

%%% Define input signal and settings

load hilbertfilter;
input = wavread('1min.wav');
fs = 8000;                              % Audio sampling rate
Ts = 1/fs;                              % Audio sampling period
winlength = 3;                          % Window size in seconds
winstart = 20;                          % Location (in seconds) to begin
hilorder = 64;                          % Order of Hilbert FIR filter
hildelay = 32;                          % Hilbert filter delay/latency
bpforder = 64;                          % Order of Bandpass Filter
bpfdelay = 32;                          % BPF delay/latency
shift = 100;                            % Desired frequency shift (in Hz)
fftorder = 500;                         % FFT length for signal analysis

%%% Prepare and analyze audio input

window = input(winstart*fs+1:winstart*fs+winlength*fs);
len = length(window);
X1 = fft(window, fftorder);
spec1 = abs(fftshift(X1)).^2';
spec1 = spec1./max(spec1);
freq = linspace(-fs/2, fs/2, fftorder);
stop1 = (abs(shift)*2)/(fs/2);
stop2 = ((fs/2)-abs(shift)*2)/(fs/2);
bandpass = fir1(bpforder, [stop1 stop2]);
bpwin_ = conv(window, bandpass);
bpwin = bpwin_(bpfdelay+1:bpfdelay+len);
bpHf = abs(freqz(bandpass, 1, fftorder/2)').^2;
bpneg = zeros(1,fftorder-length(bpHf));
for n=1:length(bpHf)
bpneg(length(bpneg)-(n-1)) = bpHf(n);   % Fold LPF frequency response
end
bpspec = [bpneg bpHf];
X2 = fft(bpwin, fftorder);
spec2 = abs(fftshift(X2)).^2';
spec2 = spec2./max(spec2);
figure;
plot(freq, spec1, 'r:'); title('Spectrum of Input Audio Signal');
xlabel('Frequency, Hz'); ylabel('Relative Power');
hold on;
plot(freq, bpspec, 'g--'); plot(freq, spec2, 'b');
legend('Original Signal','Bandpass Frequency Reponse','Filtered Signal');
hold off;

%%% SSB modulation using Hilbert filter

mhat = conv(bpwin,hilxformer);
Mr = bpwin';
Mi = mhat(hildelay+1:hildelay+len)';

% Create analytic signal and plot DFT

ssb = Mr + j*Mi;
X3 = fft(ssb, fftorder);
spec3 = abs(fftshift(X3)).^2';
spec3 = spec3./max(spec3);
hilHf = abs(freqz(hilxformer, 1, fftorder/2)').^2;
hilspec = [zeros(1, fftorder-length(hilHf)) hilHf];
figure;
plot(freq, spec2, 'r:');
title('Spectrum of Upper Sideband');
xlabel('Frequency, Hz');
ylabel('Relative Power');
hold on;
plot(freq, hilspec, 'g--');
plot(freq, spec3, 'b');
legend('Filtered Signal','Hilbert Frequency Reponse','Upper Sideband');
hold off;

%%% Modulate with "local oscillators"

n = 1:len;
sign = shift/abs(shift);
localosc = cos(2*pi*abs(shift)*n*Ts);
localosc90 = sin(2*pi*abs(shift)*n*Ts);
freqshift = localosc.*Mr - (sign)*localosc90.*Mi;
X4 = fft(freqshift, fftorder);
spec4 = abs(fftshift(X4)).^2';
spec4 = spec4./max(spec4);
figure;
plot(freq, spec2, 'k:');
title('Spectrum of Frequency Shifted Signal');
xlabel('Frequency, Hz'); ylabel('Relative Power');
hold on;
plot(freq, spec4, 'b');
legend('Filtered Signal','Shifted Signal');
hold off;

%%% To listen to your signals (at reduced volume):
% sound(window./2, fs); % orignal signal
% sound(bpwin./2, fs); % filtered signal
% sound(freqshift./2, fs); % frequency shifted signal