% class1 - script used in first ECE 201 class

%% DJIA data
% begin by importing DJIA data from file djia.csv; use "Import Data ..." in
% "File" menu

% define length 50 running average filters to be used for smoothing
M = 50;
b = ones(1,M)/M;

M2 = 15;
b2 = ones(1,M2)/M2;

% extract "Close" data ...
fid = fopen('djia_2010.csv');
C = textscan(fid,'%s %f %f %f %f %d %f', 'Delimiter',',', 'Headerlines',1);
fclose(fid);

djia  = flipud(C{5});
dates = datenum(flipud(C{1}));

... then plot
plot(dates,djia)
datetick('x','yyyy')
xlabel('Time');
ylabel('DJIA')
axis('tight')
grid

pause

% filter data for smoother signals
smooth  = conv(b,djia);
smooth2 = conv(b2,djia);

hold on;
plot(dates(51:end,1), smooth(51:length(djia)), 'r');
plot(dates(51:end,1), smooth2(51:length(djia)),'g');
hold off;

% prepare for next example
disp('press return for next example');
pause
close(gcf)

%% Music
% load sampled signal from wav file
[y,fs] = audioread('music.wav');
% display sample rate and resolution of wav file ...
fprintf('Fs: %5.0f\n',fs);
% ... then play through speaker, and ...
soundsc(y,fs)

disp('press return for plot')
pause

% plot and play again after pause
t  = (1:length(y))/fs;  % time axis
plot(t,y)
xlabel('Time (seconds)')
ylabel('Signal Amplitude')
grid on
axis('tight')

disp('press return for sound')
pause

soundsc(y,fs)

disp('press return for next example');
pause
close(gcf)

%% Heart beat signal
% read signal from wav file
[y,fs] = audioread('heart.wav');
% display sample rate and resolution
fprintf('Fs: %5.0f\n',fs);
% play through speaker
soundsc(y,fs)

disp('press return for plot')
pause

% plot and play again after pause
t = (1:length(y))/fs;    % time axis
plot(t,y)
xlabel('Time (seconds)')
ylabel('Signal Amplitude')
grid on
axis('tight')

disp('press return for sound')

pause

soundsc(y,fs)

disp('press return for next example');
pause
close(gcf)

%% Telephone dialing
fs = 8092;
tonedur = 0.25;
pausedur = 0.05;

LBPF = 131; %129;

SNR = 30; % signal-to noise ratio in dB
sigma = sqrt(0.5/10^(SNR/10));

% generate signal

xx=dtmfdial('7039931569', fs, tonedur, pausedur);
xx = xx + sigma*randn(size(xx));
tt = (0:length(xx)-1)/fs;
soundsc(xx,fs)
figure(1)
subplot(1,1,1)
plot(tt,xx)
xlabel('Time (s)')
ylabel('Signal Amplitude')
grid on
axis('tight')

disp('press return for sound')

pause

soundsc(xx,fs)

disp('press return to decode')

pause


% design BP filters
DTMFFreqs = [ 691 770 852 941 1209 1336 1477 ];

hh = zeros(length(DTMFFreqs), LBPF);

for kk=1:length(DTMFFreqs)
    hh(kk,:) = makeBPF( DTMFFreqs(kk)/fs, LBPF );
end

% filter the signal by each of the seven filters
yy = zeros(length(DTMFFreqs), length(xx)+length(hh(1,:))-1);

for kk=1:length(DTMFFreqs)
    yy(kk,:) = conv( hh(kk,:), xx );
end

% decision logic
% decompose samples into periods for each number
Nnumbers = floor( length(xx)/(fs*(tonedur+pausedur)) );
NTonePlusPause = round(fs*(tonedur+pausedur));
NPause = round(fs*pausedur);
    
% score for each tone period: sum of squares in period
score = zeros( Nnumbers, length(DTMFFreqs) );
for nn=1:Nnumbers
    Startnn = (nn-1)*NTonePlusPause + 1 + floor(LBPF/2);
    Endnn = nn*NTonePlusPause - NPause - floor(LBPF/2);
    for kk = 1:length(DTMFFreqs)
        score(nn,kk) = sum( yy(kk,Startnn:Endnn).^2 );
    end
end

% Decisions
% in each row of the score matrix find the biggest entry among the first
% four and final three columns
imax_low = zeros(1, Nnumbers );
imax_high = zeros(1, Nnumbers );

for nn=1:Nnumbers
    [ ~, imax_low(nn)] = max( score(nn, 1:4) );
    [ ~, imax_high(nn)] = max( score(nn, 5:7) );
end

% decode
% lookup table to translate numbers string into numbers
Digits = double('123456789*0#');      % table of ASCII codes for dial-pad

numbers = (imax_low-1)*3 + imax_high; % between 1 and 12
rawnumbers = Digits( numbers );       % ASCII codes for dialed numbers
decodedNumber = char( rawnumbers );   % convert to string
display( sprintf('\nDecoded Number: %s\n', decodedNumber ) );


