clear
close all
clc

load(['Mach0p8Q170.mat'])  % Load .mat file of interest
aid=5;  % Index of angle of interest from vector 'angles'
Pid=[12 13 14 15 19 29 34]; % Index of pressure taps of interest (1-35)
% FRF/PSD settings
nfft_psd=2^11;        % Number of frequncies specified
overlap_fac=0.667;  % Window overlap factor
N_overlap=floor(overlap_fac * nfft_psd);
% uses hamming window
% replace hamming(nfft_psd) with rectwin(nfft_psd) for rectangular window or
  
for j=1:numel(angles)
    % Remove Previous frequency data
    Data{j}=rmfield(Data{j},{'magPSD','freqPSD','redfreqPSD',...
        'strouhalPSD','magFFT','freqFFT','redfreqFFT','strouhalFFT'});
    for i=1:35
%         Data{j}.cPmax(i)=max(Data{j}.cP(:,i));
%         Data{j}.cPmin(i)=min(Data{j}.cP(:,i));
%         Data{j}.cPmean(i)=mean(Data{j}.cP(:,i));
        [Data{j}.magPSD(:,i) Data{j}.freqPSD(:,i)]=pwelch(Data{j}.cP(:,i)-...
            mean(Data{j}.cP(:,i)),hamming(nfft_psd),N_overlap,nfft_psd,...
            Data{j}.samp);
        
        Data{j}.redfreqPSD(:,i)=(Data{j}.freqPSD(:,i)*2*pi*(16/12))./(2*...
            Data{j}.V);
        Data{j}.strouhalPSD(:,i)=(Data{j}.freqPSD(:,i)*(16/12))./(...
            Data{j}.V);
    end
end

for j=1:numel(angles)  % FFT of full time sample (Bin for zero Hz freq excluded, way too high?)
    for i=1:35
        Data{j}.fft(:,i)=fft((Data{j}.cP(:,i)-mean(Data{j}.cP(:,i)))*...
            (Data{j}.Q/144),length(Data{j}.cP(:,i)));
        Data{j}.magFFT(:,i)=abs(Data{j}.fft(1:length(Data{j}.cP(:,i))/2+...
            1,i))*2/length(Data{j}.cP(:,i));
        Data{j}.freqFFT(:,i)=Data{j}.samp*(0:length(Data{j}.cP(:,i))/2)/...
            length(Data{j}.cP(:,i));
        Data{j}.redfreqFFT(:,i)=(Data{j}.freqFFT(:,i)*2*pi*(16/12))./(2*...
            Data{j}.V);
        Data{j}.strouhalFFT(:,i)=(Data{j}.freqFFT(:,i)*(16/12))./(...
            Data{j}.V);
    end
end

% Upper and lower x/c values corresponding to pressure taps 1-35 if needed
xocupper=[0 0.009 0.023 0.049 0.099 0.149 0.198 0.249 0.298 0.348...
    0.398 0.448 0.498 0.542 0.598 0.648 0.749 0.799 0.849 0.899 0.95 1];
xoclower=[0.012 0.027 0.103 0.203 0.303 0.403 0.503 0.602 0.652 0.702...
    0.752 0.851 0.901];
% Plot all PSDs
for i=1:length(Pid)
    figure
    loglog(Data{aid}.freqPSD(:,Pid(i)),Data{aid}.magPSD(:,Pid(i)),'r-o')
    grid on
    title(['M=' num2str(round(Data{aid}.Mach,2)) ', AoA=' num2str(angles(aid)) ', Q=' num2str(round(Data{5}.Q),-1) 'psf, Pt ' num2str(Pid(i))])
    xlabel('Frequency [Hz]')
    ylabel('Magnitude [PSI^2/Hz]')
end

%Plot all cPs
for i=1:length(Pid)
    figure
    plot(Data{aid}.time,Data{aid}.cP(:,Pid(i)),'b-')
    grid on
    title(['M=' num2str(round(Data{aid}.Mach,2)) ', AoA=' num2str(angles(aid)) ', Q=' num2str(round(Data{5}.Q),-1) 'psf, Pt ' num2str(Pid(i))])
    xlabel('Time [seconds]')
    ylabel('Cp')
end