%=======================================================================================================================
% THDN_TOP.m
%         
% Demonstrator of the function THDN
%
%
% THDN computes the total harmonic distortion + noise ratio ratio of a SINEWAVE 
%  signal corrupted with noise from different sources (jitter, saturation, pink noise, ..)
%
% Program validated with two versions of Matlab:
%  Matlab Version 7.2 (R2006a), Signal Processing Toolbox Version 6.5 (R2006a)
%  Matlab Version R2015
%  
close all; clear all;

select_example = 3;
%
% EXPECTED RESULTS FROM THE THREE EXAMPLES
%
% EXAMPLE1
% THDN     =  -67.718  =   11.249 bits
% THDN_CCIR=  -68.225  =   11.333 bits
%
% EXAMPLE2
% THDN     = -124.436  =   20.670 bits
% THDN_Aw  = -132.327  =   21.981 bits
%
% EXAMPLE3 (Calibration)
% ORIGINAL = -144.500  =   24.003 bits
% THDN     = -144.425  =   23.991 bits
% THDN_Aw  = -147.104  =   24.436 bits



%======= EXAMPLE OF VOICE FILE ANALYSIS =======
if select_example == 1
    % THD of PDM data
    [data_in, Fs, NBITS] = wavread('VOICEBAND.WAV');  
    %[data_in, Fs] = audioread('VOICEBAND.WAV');   % instruction to be used starting from Matlab2015
    Selection_A0_CCIR1 = 1;     
    [thd, thd_weighted, ratio_to_full_band] = THDN (data_in, Fs, 300, 3300, Selection_A0_CCIR1);
end

%======= EXAMPLE OF PDM FILE ANALYSIS =======
if select_example == 2
    %
    % THD of PDM data recorded at 6.144MHz in a text file
    %  Estimate the THD+N with and without A-weigthed law
    %
    [fid_1] = fopen('PDM6144.TXT','r');
    [pdm_data] = fscanf (fid_1,'%d'); fclose(fid_1);
    Selection_A0_CCIR1 = 0;     
    [thd, thd_weighted, ratio_to_full_band] = THDN (pdm_data, 6144000, 20, 20000, Selection_A0_CCIR1);
    
    pdm_modulation_index = ratio_to_full_band                           % -13.436 dB in the example
end

%======= CHECK AT MATLAB ACCURACY BOUNDARIES =========
if select_example == 3
    NBITS = 24                                                          % 1bit .. 26bits audio accuracy
    data_length = 500000;                                               % take enough samples for this measurement
    Fs = 48000;                                                         % artificial data at 48kHz
    t1 = 2^NBITS * sqrt(2) * sin(2*pi*[1:data_length]/(Fs/997))';       % sinewave signal, 24bits 
    t2 = (2^0 *(2/0.57736)*(0.5- rand(1,data_length)));                 % noise
    
    EX1 = sqrt(t1'*t1) / sqrt(data_length);                             % signal energy
    EX2 = sqrt(t2*t2') / sqrt(data_length);                             % noise energy 
    fprintf('ORIGINAL = %8.3f  = %8.3f bits\n', 20*log10(EX2/EX1),  abs(20*log10(EX1/EX2))/6.02);
    data_in = t1'+t2; 
    Selection_A0_CCIR1 = 0;
    [thd, thd_weighted, ratio_to_full_band] = THDN (data_in, Fs, 1, 24000, Selection_A0_CCIR1);    
end


%
%======= PRINT RESULTS =======
%
fprintf('THDN     = %8.3f  = %8.3f bits\n', thd,  abs(thd)/6.02);
if Selection_A0_CCIR1==0
    fprintf('THDN_Aw  = %8.3f  = %8.3f bits\n', thd_weighted,  abs(thd_weighted)/6.02);
else
    fprintf('THDN_CCIR= %8.3f  = %8.3f bits\n', thd_weighted,  abs(thd_weighted)/6.02);
end

