% In this program, I have a single observation (y_1, y_2) from a bivariate normal 
% distribution with known covariance matrix and unknown mean vector 
% (th(1), th(2)). The covariance matrix has ones on its principal 
% diagonal and rho elsewhere.  
% Author : ChangHyun Kim 
% Date : Feb 1st 2007
%
clear; close; clc;
y_1 = 0; y_2 = 0;  % the observed values of y_1 and y_2
rho = .8;          % correlation coefficient 
M = 120;           % number of samples to be drawn
gibbssample = zeros(M,2);   % create an Mx2 matrix to store draws
gibbssample(1,1)=2.0; gibbssample(1,2)=-2.0;  % initial values of parameters
% Now we begin gibbssample sampling, using formulas based on gibbssamplee 
% bivariate normal distribution.
for ii=2:M                   
    
 gibbssample(ii,1)=y_1+rho*(gibbssample(ii-1,2)-y_2)+sqrt(1-rho^2)*randn(1,1);
 gibbssample(ii,2)=y_2+rho*(gibbssample(ii,1)-y_1)+sqrt(1-rho^2)*randn(1,1);
end;
% Plots
% plot(gibbssample(:,1),gibbssample(:,2),'-')   % plot pagibbssample gibbssamplerough parameter space
% xlabel('\gibbssampleeta_1','fontsize',14)
% ylabel('\gibbssampleeta_2','fontsize',14)
% legend('gibbs sample 1D two arrays');
% pause
scatter(gibbssample(:,1),gibbssample(:,2),'filled') % plot whole samples
xlabel('\gibbssampleeta_1','fontsize',14)
ylabel('\gibbssampleeta_2','fontsize',14)
legend('gibbs sample 1D two arrays');

% Sound Alogorithm for A major tonal music

fs = 44100;    % sampling frequency
t = 1:1:22050; % each (event)note duration length 0.5sec
Key = [1 2^(2/12) 2^(4/12) 2^(5/12) 2^(7/12) 2^(9/12) 2^(11/12) 2]; % A major scale chord
numNotes = 120; % number of notes (for 1 min piece, the same as number of samples)
y = zeros(1,fs/2); 
yy = zeros(1,fs*60);
% 
mapping2 = 0;
for j=1:numNotes,
    if gibbssample(j,1) < -2.0 | gibbssample(j,1) >= 2
        y = abs(gibbssample(j,2)).*randn(1,fs/2);
        if mapping2 == 1
            y = abs(gibbssample(j,2)).*zeros(1,fs/2);
        end
    elseif gibbssample(j,1) >= -2.0 & gibbssample(j,1) < -1.5
            freq = 440*Key(1);
            if mapping2 == 1
                freq = 440*Key(4);
            end
            y = abs(gibbssample(j,2)).*cos(2*pi*freq/fs*t);
    elseif gibbssample(j,1) >= -1.5 & gibbssample(j,1) < -1.0
            freq = 440*Key(2);
            y = abs(gibbssample(j,2)).*cos(2*pi*freq/fs*t);
    elseif gibbssample(j,1) >= -1.0 & gibbssample(j,1) < -0.5
            freq = 440*Key(3);
            y = abs(gibbssample(j,2)).*cos(2*pi*freq/fs*t);
    elseif gibbssample(j,1) >= -0.5 & gibbssample(j,1) < 0.0
            freq = 440*Key(4);
            if mapping2 == 1
                freq = 440*Key(1);
            end
            y = abs(gibbssample(j,2)).*cos(2*pi*freq/fs*t);
    elseif gibbssample(j,1) >= 0.0 & gibbssample(j,1) < 0.5
            freq = 440*Key(5);
            if mapping2 == 1
                freq = 440*Key(8);
            end
            y = abs(gibbssample(j,2)).*cos(2*pi*freq/fs*t);
    elseif gibbssample(j,1) >= 0.5 & gibbssample(j,1) < 1.0
            freq = 440*Key(6);
            y = abs(gibbssample(j,2)).*cos(2*pi*freq/fs*t);
    elseif gibbssample(j,1) >= 1.0 & gibbssample(j,1) < 1.5
            freq = 440*Key(7);
            y = abs(gibbssample(j,2)).*cos(2*pi*freq/fs*t);
    elseif gibbssample(j,1) >= 1.5 & gibbssample(j,1) < 2.0
            freq = 440*Key(8);
            if mapping2 == 1
                freq = 440*Key(5);
            end
            y = abs(gibbssample(j,2)).*cos(2*pi*freq/fs*t);
    end
    if j == 1
        yy = y;
    end
    if j ~= 1
        yy = [yy y];
    end
end
%sound(yy, fs);
wavwrite(yy,fs,8,'gibbssampling.wav');