
%takes a filename to output the error to
%        an array of data (one dimensional)
%        an array of dates 
%pmoption - plot models option, 'p' to plot
%outputs array of models same length as the original data array

function [out,selthresh]=DCTthresholdSelection(A,numtrpoints,n,plotoption)  

data = A(1:numtrpoints);
Cm =[];
clear i;
B = dct(data);
[sorted_data,ind] = sort(abs(B));   %store DCT coefficients in order of increasing magnitude
RMST = [];                          %training RMS (RMS(model-data))
RMSPD = [];                         %prediction RMS vs Data
RMSPM = [];                         %prediction RMS vs model
SUMRMS = [];
ld=length(sorted_data);

%for i=ld-19:-n:ld-58
for i=ld:-n:1
    B = dct(data);
    temp = find(abs(B)<sorted_data(i));  %find all coefficients that are smaller than given coeff (remove smaller ones first) 
    B(temp) = zeros(size(temp));
    new_model = idct(B);
    predlistmodel = predstd(new_model,28,n,[]);    
    
    %RMS between model and data
    rmstr = norm(new_model(1:numtrpoints)-data(1:numtrpoints))/sqrt(numtrpoints);
    RMST = [RMST, rmstr];
    
    % RMS between data and prediction
    rmspd = norm(data(29:numtrpoints)-predlistmodel(:,2))/sqrt(numtrpoints);
    RMSPD = [RMSPD, rmspd];
    
    % RMS between prediction and the model
    rmspm = norm(predlistmodel(:,1)-predlistmodel(:,2))/sqrt(numtrpoints);
    RMSPM = [RMSPM, rmspm];
    
    SUMRMS = [SUMRMS, rmstr+rmspd+rmspm];
end;
indmin = find(SUMRMS == min(SUMRMS));

B = dct(data);
temp = find(abs(B)<sorted_data(ld-18-indmin));  %find all coefficients that are smaller than given coeff (remove smaller ones first) 
B(temp) = zeros(size(temp));
out = idct(B);

fid = fopen('DCTcoughMINSum.rms','a+');
fprintf(fid,'%d %f %f\n',indmin,sorted_data(ld-18-indmin),SUMRMS(indmin));
fclose(fid);

%selthresh = sorted_data(ld-18-indmin);
indmin
selthresh = sorted_data(ld-indmin)

if (plotoption == 'p')
    figure;
%    subplot(3,1,1);
    plot(RMST,'--','Color','black');
    hold on;
    plot(RMSPD);
    plot(RMSPM,'-.','Color','black');
    legend('RMS(model,data)','RMS(pred,data)','RMS(pred,model)');
    xlabel('all coefficients smaller than the coefficient with this index number are removed');
    ylabel('rms error');
    title(['Predicting ' num2str(n) ' values in the future from ' num2str(numtrpoints) ' observations ']);
 %   fprintf('Prediction RMS is %f with threshold %f.  Removed coefficients %d\n',minp,threshold,coefrem);
 %   fprintf('Positive prediction RMS is %f with threshold %f.  Removed coefficients %d\n',minpp,thresholdp,coefremp);
    hold off;
%    subplot(3,1,2);
   figure
   plot(SUMRMS);
   legend('sum(RMSs)');
end;
