%% This routine calculates the spectral pitch similarity of any given pitch 
% class set with a series of 'probe' pitch classes or chords. The spectral
% pitch vectors parameterized by rolloff (rho) and smoothing width (sigma).

% Author: Andrew J. Milne
% Version: 1.0

clear all

global virtual

%% Options
q = 1200;
period = 1200;
metric = 'cosine';
p = 1;
sigma = 5.9;
win_len = 12;
virtual = 0;
abs_or_rel = 'abs';
r = 1;
exp_or_sal = 'exp'

num_harmonics = 12;
rolloff = 0.67;
remove_dupes = 0;
normed_corr = 1;

%%
y_Maj = [6.35 2.23 3.48 2.33 4.38 4.09 2.52 5.19 2.39 3.66 2.29 2.88]'; % p.30 - these are the major hierarchies
y_min = [6.33 2.68 3.52 5.38 2.60 3.53 2.54 4.75 3.98 2.69 3.34 3.17]'; % p.30 - these are the minor hierarchies
y = [y_Maj y_min]-1;

%pc_set = [5 9 0 4 7 11];
% pc_set = [5 9 0 7 11 2];
%pc_set = [0 4 7];
%pc_set = [2 5 9  4 8 11  9 0 4] % d - E - a
%pc_set = [7 11 2  5 8 0  0 3 7] % G - f - c
% pc_set = [5 9 0  7 11 2] % F - G  
% pc_set = [5 9 0  7 10 2] % F - g  
% pc_set = [2 5 9  7 11 2] % d - G  
% pc_set = [2 5 9  7 10 ] % d - g  
% pc_set = [2 5 9  7 11 2  0 4 7] % d - G - C
% pc_set = [2 5 9  7 10 2  0 4 7] % d - g - C
% pc_set = [1 5 8  7 11 2] % Db - G
% pc_set = [1 5 8  7 11 2  0 3 7] + 1 % - Db - G - c
% pc_set = [10 2 5  11 2 6] % Db - G
% pc_set = [3 7 10  7 11 2] % Eb - G
% pc_set = [11 2 5 0 4 7] % - b dim + C
% pc_set = [11 2 5 7 11 2] % - b dim + G

% pc_set = [0 4 7] % major triad
% pc_set = [0 3 7] % major triad

% pc_set = [0 2 4 7 9] % pentatonic hexachord
% pc_set = [0 2 4 5 7 9] % diatonic hexachord
pc_set = [0 2 4 5 7 9 11] % diatonic
pc_set = [0 0 0 9 11 11] % maj scale segment
pc_set = [0 3] % min scale segment
% pc_set = [0 2 3 5 7 8 11] % - harm minor
% pc_set = [0 2 4 5 7 8 11] % - harm major
% pc_set = [0 2 3 5 7 9 11] % - mel minor
% pc_set = [0 2 3 6 7 8 11] % - dble harm minor
% pc_set = [0 2 4 6 8 10] % - whole tone
% pc_set = [0 2 3 5 6 8 9 11] % - diminished
% pc_set = [0 1 3 5 7 9 11] % - whole plus one

% pc_set = [7 11 2 5] % G7
% pc_set = [0 2 4 5 7 7 11]; % C + G7
%pc_set = [0 2 4 5 7 11]; % C + Bdim
%pc_set = [0 2 4 7 11]; % C + G
% pc_set = [0 4 5 7 11 2];
% pc_set = [0 2 4 5 7 9]
% pc_set = [0 4 7 ([0 4 7] + 4)]
% IS IT POSSIBLE TO MAKE A SCALE WITHOUT TRITONES WHERE CHORDS HAVE
% HETERGENOUS FITS? OF COURSE, A TRITONE DOES NOT GUARANTEE HETEROGENAIETY,
% BUT IS IT A NECESSARY CONDITION? IT IS NOT A NECESSARY CINDITION - E.G., 
% [0 4 7 8 11]
% pc_set = [0 2 4 5 7 8 9 11]; % flamenco


% pc_set = [0 2 4 5 7 9 11  4 ] % diatonic + e
% pc_set = [0 2 4 5 7 9 11  2] % diatonic + G


% pc_set = [5 9 0  7 11 2  0 4 7] % F - G - C
% pc_set = [5 9 0  4 7 11  0 4 7] % F - e - C
% pc_set = [5 9 0  4 7 11] % F - e  
% pc_set = [5 9 0  4 7 11  0 4 7] % F - e - C  

%% Correlations between pc_set and the 12 major and all 12 minor tonal hierarchies
pc_set = mod(pc_set,12);
if remove_dupes==1
    pc_set = unique(pc_set);
end
pc_set

% Embed pc_set into an indicator vector
pc_set_ind = zeros(1,12);
pc_set_ind(pc_set+1) = 1;

if normed_corr==0
    for j=1:2
        for i=1:12
            keycentre(j,i) ...
                = pc_set_ind*circshift(y(:,j),i-1);
        end
    end
    keycentre; % This shows the "unnormalized correlation" of each tonal centre
    % with the pitch class set. Major tonalities in the top row, minor in the
    % bottom row
end

if normed_corr==1
    for j=1:2
        for i=1:12
            keycentre(j,i) ...
                = corr(pc_set_ind', circshift(y(:,j),i-1));
        end
    end
    keycentre; % This shows the correlation of each tonal centre
    % with the pitch class set. Major tonalities in the top row, minor in the
    % bottom row
end

%% plots
figure(203),subplot(1,2,1)
bar(keycentre(1,:),'FaceColor',[0.2 0.2 0.2])
if normed_corr==1
    axis([0,13,-1,1])
elseif normed_corr==0
    axis([0,13,0.9*min(min(keycentre)),1.1*max(max(keycentre))])
end
set(gca,'XTickLabel',0:12)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])

figure(203),subplot(1,2,2)
bar(keycentre(2,:),'FaceColor',[0.2 0.2 0.2])
if normed_corr==1
    axis([0,13,-1,1])
elseif normed_corr==0
    axis([0,13,0.9*min(min(keycentre)),1.1*max(max(keycentre))])
end
set(gca,'XTickLabel',0:12)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])

%% options 
% for n-TETs
low_limit = 0;
high_limit = 1200;
inc = 100;
shifts = low_limit:inc:high_limit-inc;

len_shifts = length(shifts);
Dist = zeros(1,len_shifts);

%% Generate reference pc set's spectrum
SpectrumX = zeros(1, num_harmonics);
for j = 1:num_harmonics
    SpectrumX(j) = 1200*log2(j); % !!!!!!!
end

SpectrumX_w = zeros(1, num_harmonics);
for j = 1:num_harmonics
    SpectrumX_w(j) = 1/(j^rolloff);
end

for i=1:length(pc_set)
    SpectrumX_pc_set(:,i) = (SpectrumX + 100*pc_set(i))';
end
SpectrumX_pc_set = SpectrumX_pc_set(:);
SpectrumX_w_pc_set = repmat(SpectrumX_w,[1 length(pc_set)]);

SpectrumY_base = SpectrumX;
SpectrumY_w = SpectrumX_w;

%% Generate continuum tones and send to embeddings
for i = 1:len_shifts
    shift_by = shifts(i);
    SpectrumY = SpectrumY_base+shift_by;
    spec_dist(i) ...
        = ExpectationMetrics(period,q,metric,p,sigma,win_len,...
        SpectrumX_pc_set,SpectrumX_w_pc_set,SpectrumY,SpectrumY_w,'abs',1,exp_or_sal);
end
spec_sim = 1 - spec_dist % this gives the spectral pitch similarity of all 12 
% pitch classes and pc_set.

tones_in = zeros(1,12);
tones_in(pc_set+1) = 1;
spec_sim_in = tones_in.*spec_sim % spec_sim of tones in the scale
spec_sim_out = (1-tones_in).*spec_sim % spec_sim of tones not in the scale

figure(209)
bar(spec_sim_out,'FaceColor',[0.6 0.6 0.6])
hold on
bar(spec_sim_in,'FaceColor',[0.2 0.2 0.2])
hold off
axis([0,13,0,1])
set(gca,'XTickLabel',0:12)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])

all_played = mean(spec_sim(pc_set+1)) % this gives the mean spectral pitch 
% similarity of each tone with the tones in pc_set

%% This calculates the fit of triads to pc_set directly (by comparing their 
% spectra) The two methods produce very similar results, the overall values
% of this method are a bit higher but proportional compared to the previous

SpectrumY_base = [SpectrumX 400+SpectrumX 700+SpectrumX; ...
                 SpectrumX 300+SpectrumX 700+SpectrumX];
SpectrumY_w = repmat(SpectrumX_w,[1 3]);
for j=1:2
    for i = 1:len_shifts
        shift_by = shifts(i);
        SpectrumY = SpectrumY_base(j,:)+shift_by; %%%%% !!!!!!!
        spec_dist(j,i) ...
            = ExpectationMetrics(period,q,metric,p,sigma,win_len,...
            SpectrumX_pc_set,SpectrumX_w_pc_set,SpectrumY,SpectrumY_w,'abs',1,exp_or_sal);
    end
end
spec_sim = 1 - spec_dist % this gives the spectral pitch similarity of all 12 
% pitch classes and pc_set.

% This gives the indices of the roots of the triads that are major (first row)
% or minor (second row)
chord_in = zeros(2,12);
for j=1:2
    for i=1:length(pc_set)
        if sum(ismember([0 5-j 7], mod(pc_set-pc_set(i),12))) == 3
            chord_in(j,pc_set(i)+1) = 1;
        end
    end
end
            

chord_fit_in = spec_sim.*chord_in % how well does every major minor triad 
% that is in the scale fit that scale
chord_fit_out = spec_sim.*(1-chord_in) % how well does every major minor triad 
% that is not in the scale fit that scale

figure(212)
hold off
%bar(chord_fit_out(1,:),'FaceColor',[0.6 0.6 0.6]);figure(gcf)
%bar(chord_fit_out(2,:),'FaceColor',[0.6 0.6 0.6]);figure(gcf)
bar(chord_fit_in(2,:),'FaceColor',[0.6 0.6 0.6]);figure(gcf)
hold on
bar(chord_fit_in(1,:),'FaceColor',[0.2 0.2 0.2]);figure(gcf)
axis([0,13,0,1])
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])
set(gca,'XTickLabel',0:12)

%% Non-12-TET tunings
clear all

%% Options
q = 1200;
period = 1200;
metric = 'cosine';
p = 1;
sigma = 5.9;
win_len = 12;
virtual = 0;
abs_or_rel = 'abs';
r = 1;
exp_or_sal = 'exp'

num_harmonics = 12;
rolloff = 0.67;
remove_dupes = 0;
normed_corr = 1;

%% Uncomment whichever non-12-TET tunings/scales you want to test ...

%% Pythagorean (note here I am tetsing the fit of perfect fifth dyads, since these would have been common finals - not major or minor triads)
%{
pc_set = [0 204 408 498 702 906 ] % pyth major hexachord

maj_triad = [0 702]; % pyth dyad
min_triad = [0 702];

qc_fifth = 1200*log2(3/2)
shifts = round([0 -5 7 2 -3 9 4 -1 -6 6 1 -4 8 3 -2 10 5] * qc_fifth); % Pyth 17
shifts = mod(shifts,1200) % this is the "chromatic" scale

    in_scale = [1 0 0 1 0 0 1 1 0 0 1 0 0 1 0 0 0 ]; % diatonic major

maj_in_scale = [1 0 0 1 0 0 0 1 0 0 1 0 0 1 0 0 0 ]; % maj chord roots in the scale
min_in_scale = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]; % min chord roots in the scale

maj_triad = [0 408 702]; % pyth maj triad
min_triad = [0 294 702]; % pyth min triad

maj_in_scale = [1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 ]; % maj chord roots in the scale
min_in_scale = [0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 ]; % min chord roots in the scale
%}

%% 1/4-comma meantone
%{
pc_set = [0 193 386 503 697 890 1083] % meantone diatonic major scale

maj_triad = [0 386 697]; % meantone major triad
min_triad = [0 310 697]; % meantone minor triad

qc_fifth = 1200*log2((3/2)/((81/80)^(1/4)))
shifts = round([0 7 -5 2 9 -3 4 11 -1 6 -6 1 8 -4 3 10 -2  5 -7] * qc_fifth); % meantone 19 1/4-comma
shifts = mod(shifts,1200) % this is the "chromatic" scale

in_scale = [1 0 0 1 0 0 1 0 1 0 0 1 0 0 1 0 0 1 0 ]; % diatonic major
% in_scale = [1 0 0 1 0 1 0 0 1 0 0 1 0 1 0 0 0 1 0 ]; % harmonic minor

maj_in_scale = [1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 ]; % maj chord roots in the scale
min_in_scale = [0 0 0 1 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 ]; % min chord roots in the scale
%}

%% porcupine in 22-tet
%{
maj_triad = [0 382 709]; % porcupine major triad
min_triad = [0 327 709]; % porcupine minor triad

eq22 = 1200/22
shifts = round([0:21] * eq22); % 
shifts = mod(shifts,1200) % this is the "chromatic" scale

in_scale = [1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 ]; % 1L, 6s
not_in_scale = 1-in_scale;
pc_set = shifts;
pc_set(logical(not_in_scale)) = []; % porcupine 1L, 6s scale

maj_in_scale = [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 ];
min_in_scale = [0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 ];
%}

%% srutal in 22-tet
%{
maj_triad = [0 382 709]; % srutal major triad
min_triad = [0 327 709]; % srutal minor triad

eq22 = 1200/22
shifts = round([0:21] * eq22); % 
shifts = mod(shifts,1200) % this is the "chromatic" scale

in_scale = [1 0 1 0 1 0 1 0 1 0 0 1 0 1 0 1 0 1 0 1 0 0 ]; % 1L, 6s
not_in_scale = 1-in_scale;
pc_set = shifts;
pc_set(logical(not_in_scale)) = []; % srutal 1L, 4s scale

maj_in_scale = [0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 1 0 1 0 1 0 0 ];
min_in_scale = [1 0 1 0 1 0 0 0 0 0 0 1 0 1 0 1 0 0 0 0 0 0 ];
%}

%% magic in 22-tet
%{
maj_triad = [0 382 709]; % srutal major triad
min_triad = [0 327 709]; % srutal minor triad

eq22 = 1200/22
shifts = round([0:21] * eq22); % 
shifts = mod(shifts,1200) % this is the "chromatic" scale

in_scale = [1 1 1 0 0 0 0 1 1 1 0 0 0 0 1 1 1 1 0 0 0 0 ]; % 3L, 7s
not_in_scale = 1-in_scale;
pc_set = shifts;
pc_set(logical(not_in_scale)) = [];

maj_in_scale = [0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 ];
min_in_scale = [0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 ];
%}

%% Uncomment to graph results for microtonal scales
%{
len_shifts = length(shifts);
Dist = zeros(1,len_shifts);

%% Generate reference pc set's spectrum
SpectrumX = zeros(1, num_harmonics);
for j = 1:num_harmonics
    SpectrumX(j) = 1200*log2(j);
end

SpectrumX_w = zeros(1, num_harmonics);
for j = 1:num_harmonics
    SpectrumX_w(j) = 1/(j^rolloff);
end

for i=1:length(pc_set)
    SpectrumX_pc_set(:,i) = (SpectrumX + pc_set(i))';
end
SpectrumX_pc_set = SpectrumX_pc_set(:);
SpectrumX_w_pc_set = repmat(SpectrumX_w,[1 length(pc_set)]);

SpectrumY_base = SpectrumX;
SpectrumY_w = SpectrumX_w;

%% Generate continuum tones and send to embeddings
spec_dist = zeros(1,1200);
for i = 1:len_shifts
    shift_by = shifts(i);
    SpectrumY = SpectrumY_base+shift_by
    spec_dist(shift_by+1) ...
        = ExpectationMetrics(period,q,metric,p,sigma,win_len,...
        SpectrumX_pc_set',SpectrumX_w_pc_set',SpectrumY,SpectrumY_w,'abs',1,exp_or_sal);
end
spec_sim = 1 - spec_dist; % this gives the spectral pitch similarity of all 12 
% pitch classes and pc_set.

% tones_in = zeros(1,1200);
% tones_in(pc_set+1) = 1;
% spec_sim_in = tones_in.*spec_sim % spec_sim of tones in the scale
% spec_sim_out = (1-tones_in).*spec_sim % spec_sim of tones not in the scale

out_out = 1:1200;
out_out = setdiff(out_out,shifts+1);

spec_sim_out = spec_sim;
spec_sim_out(out_out) = []

spec_sim_in = spec_sim_out.*in_scale % diatonic major

figure(309)
bar(spec_sim_out,'FaceColor',[0.6 0.6 0.6])
hold on
bar(spec_sim_in,'FaceColor',[0.2 0.2 0.2])
hold off
axis([0,length(shifts)+1,0,1])
set(gca,'XTick',1:length(shifts))
set(gca,'XTickLabel',0:length(shifts)-1)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])

% theta = 2*pi*(0:100:1200)/period;
% figure(2009), mmpolar(theta', [spec_sim spec_sim(1)], 'Style', 'compass', ...
%     'RLimit', [0 1], ...
%     'RTickLabelVisible', 'off', ...
%     'TTickValue', [0 30 60 90 120 150 180 210 240 270 300 330], ...
%     'TTickLabel', {'0' '1' '2' '3' '4' '5' '6' '7' '8' '9' '10' '11'}, ...
%     'FontSize', 14)

all_played = mean(spec_sim(pc_set+1)) % this gives the mean spectral pitch 
% similarity of each tone with the tones in pc_set


%% This calculates the fit of triads to pc_set directly (by comparing their 
% spectra) The two methods produce very similar results, the overall values
% of this method are a bit higher but proportional compared to the previous

if length(maj_triad)==3 && length(min_triad)==3
    SpectrumY_base = [SpectrumX maj_triad(2)+SpectrumX maj_triad(3)+SpectrumX; ...
                     SpectrumX min_triad(2)+SpectrumX min_triad(3)+SpectrumX];
    SpectrumY_w = repmat(SpectrumX_w,[1 3]);
elseif length(maj_triad)==2 && length(min_triad)==2
    SpectrumY_base = [SpectrumX maj_triad(2)+SpectrumX; ...
                     SpectrumX min_triad(2)+SpectrumX];
    SpectrumY_w = repmat(SpectrumX_w,[1 2]);
end
    
    
for j=1:2
    for i = 1:len_shifts
        shift_by = shifts(i);
        SpectrumY = SpectrumY_base(j,:)+shift_by;
        spec_dist(j,shift_by+1) ...
            = ExpectationMetrics(period,q,metric,p,sigma,win_len,...
            SpectrumX_pc_set,SpectrumX_w_pc_set,SpectrumY,SpectrumY_w,'abs',1,exp_or_sal);
    end
end
spec_sim = 1 - spec_dist % this gives the spectral pitch similarity of all 12 
% pitch classes and pc_set.

spec_sim_out = spec_sim;
spec_sim_out(:,out_out) = []

%spec_sim_in = spec_sim_out.*[1 0 0 1 0 0 1 0 1 0 0 1 0 0 1 0 0 1 0; 1 0 0 1 0 0 1 0 1 0 0 1 0 0 1 0 0 1 0]

chord_maj = spec_sim_out.*[maj_in_scale; zeros(1,length(maj_in_scale)) ]
chord_min = spec_sim_out.*[zeros(1,length(min_in_scale)); min_in_scale ]

for i=1:length(chord_maj)
    if chord_maj(1,i) >= chord_min(2,i)-0.0000000000001
        max_maj(i) = 1;
        max_min(i) = 0;
    else
        max_maj(i) = 0;
        max_min(i) = 1;
    end
end

chord_maj = chord_maj.*[max_maj;max_maj]
chord_min = chord_min.*[max_min;max_min]

figure(312)
hold off
bar(chord_min(2,:),'FaceColor',[0.6 0.6 0.6]);figure(gcf)
hold on
bar(chord_maj(1,:),'FaceColor',[0.2 0.2 0.2]);figure(gcf)
axis([0,length(shifts)+1,0,1])
set(gca,'XTick',1:length(shifts))
set(gca,'XTickLabel',0:length(shifts)-1)
set(gca,'PlotBoxAspectRatio',[1.618034,1,1])

%}