From edb3900153484ce521bd2c492f69c63257f71177 Mon Sep 17 00:00:00 2001 From: mike-dixon Date: Fri, 9 Aug 2024 18:03:30 -0600 Subject: [PATCH] Smoothing factor analysis --- .../timeSeries/noisePeaks_smoothingTest.m | 7 +- .../timeSeries/run_smoothingTest.m | 4 +- .../timeSeries/smoothingPlots.m | 236 ++++++++++++++++++ 3 files changed, 242 insertions(+), 5 deletions(-) create mode 100644 projDir/qc/dataProcessing/timeSeries/smoothingPlots.m diff --git a/projDir/qc/dataProcessing/timeSeries/noisePeaks_smoothingTest.m b/projDir/qc/dataProcessing/timeSeries/noisePeaks_smoothingTest.m index 16cb321e..f97d25ad 100644 --- a/projDir/qc/dataProcessing/timeSeries/noisePeaks_smoothingTest.m +++ b/projDir/qc/dataProcessing/timeSeries/noisePeaks_smoothingTest.m @@ -1,4 +1,4 @@ -function [err,resid]=noisePeaks_smoothingTest(specDB,velIn,data,widthC,err,resid,figdir,plotTime) +function [err,resid]=noisePeaks_smoothingTest(specDB,velIn,data,widthC,aircVel,err,resid,figdir,plotTime) % Find mean noise and noise threshold with following % Hildebrand and Sekhon, 1974 https://doi.org/10.1175/1520-0450(1974)013%3C0808:ODOTNL%3E2.0.CO;2 % Adjust spectra so they fit in the boundaries @@ -64,7 +64,8 @@ % VEL meanVel=sum(sigInLin.*testVel,'omitmissing')/sum(sigInLin,'omitmissing'); - filterAt=8; + %filterAt=8; + filterAt=round(0.00022396*aircVel^2-0.10542*aircVel+18.132); % Correct for aircraft width [err,errCat,sigWidthCorr,sigFiltered,signalIn1,signalIn2,sigFiltered1,sigFiltered2,inds1,inds2]= ... @@ -123,7 +124,7 @@ xlim([testVel(1),testVel(end)]); - legend('Original signal','Filtered','Filtered width corrected') + legend('Original signal','Filtered',['Filtered width corrected (',num2str(filterAt),' nz)']); grid on box on diff --git a/projDir/qc/dataProcessing/timeSeries/run_smoothingTest.m b/projDir/qc/dataProcessing/timeSeries/run_smoothingTest.m index c920d950..10a074aa 100644 --- a/projDir/qc/dataProcessing/timeSeries/run_smoothingTest.m +++ b/projDir/qc/dataProcessing/timeSeries/run_smoothingTest.m @@ -276,7 +276,7 @@ % This step removes the noise, de-aliases, (and corrects for % spectral broadening) [err,resid]=noisePeaks_smoothingTest(specPowerDB.V, ... - momentsTimeOne.velRawDeAliased(:,ii),dataThis,widthCorrDelta(cfInd),err,resid,figdir,plotTime); + momentsTimeOne.velRawDeAliased(:,ii),dataThis,widthCorrDelta(cfInd),velAircraft(cfInd),err,resid,figdir,plotTime); velAirc=cat(2,velAirc,repmat(velAircraft(cfInd),1,size(err,2)-size(velAirc,2))); end @@ -516,7 +516,7 @@ set(gcf,'PaperPositionMode','auto') print(f1,[figdir,project,'_smoothingAnalysis_everyOther_aircraftSpeed'],'-dpng','-r0'); -save([figdir,project,'_smoothingAnalysis_everyOther_aircraftSpeed.mat'],'velAircAll','errAll','numZero'); +save([figdir,project,'_smoothingAnalysis_everyOther_aircraftSpeed.mat'],'edgesHalf','peakNum'); %% Plot cases diff --git a/projDir/qc/dataProcessing/timeSeries/smoothingPlots.m b/projDir/qc/dataProcessing/timeSeries/smoothingPlots.m new file mode 100644 index 00000000..050a24d1 --- /dev/null +++ b/projDir/qc/dataProcessing/timeSeries/smoothingPlots.m @@ -0,0 +1,236 @@ +% Analyze HCR time series +clear all; +close all; + +addpath(genpath('~/git/HCR_configuration/projDir/qc/dataProcessing/')); + +figdir='/scr/virga1/rsfdata/projects/spicule/hcr/qc1/cfradial/v1.2_full/spectralMoments/smoothingTest/projectComparison/'; + +%% Speed bins +socrates=load('/scr/virga2/rsfdata/projects/socrates/hcr/qc3/cfradial/v3.2_full/spectralMoments/smoothingTest/socrates_smoothingAnalysis_everyOther_aircraftSpeed.mat'); +otrec=load('/scr/sleet2/rsfdata/projects/otrec/hcr/qc3/cfradial/v3.2_full/spectralMoments/smoothingTest/otrec_smoothingAnalysis_everyOther_aircraftSpeed.mat'); + +%% Combine +socBoth=cat(2,socrates.edgesHalf',socrates.peakNum); +socBoth(any(isnan(socBoth),2),:)=[]; +otBoth=cat(2,otrec.edgesHalf',otrec.peakNum); +otBoth(any(isnan(otBoth),2),:)=[]; + +comb=cat(1,socBoth,otBoth); +comb=sortrows(comb); + + +%% Fits + +[pLs,sLs]=polyfit(socBoth(:,1),socBoth(:,2),1); +[yLs,dLs]=polyval(pLs,socBoth(:,1),sLs); + +[pCs,sCs]=polyfit(socBoth(:,1),socBoth(:,2),2); +[yCs,dCs]=polyval(pCs,socBoth(:,1),sCs); + +[pLo,sLo]=polyfit(otBoth(:,1),otBoth(:,2),1); +[yLo,dLo]=polyval(pLo,otBoth(:,1),sLo); + +[pCo,sCo]=polyfit(otBoth(:,1),otBoth(:,2),2); +[yCo,dCo]=polyval(pCo,otBoth(:,1),sCo); + +[pLb,sLb]=polyfit(comb(:,1),comb(:,2),1); +[yLb,dLb]=polyval(pLb,comb(:,1),sLb); + +[pCb,sCb]=polyfit(comb(:,1),comb(:,2),2); +[yCb,dCb]=polyval(pCb,comb(:,1),sCb); + +%% Scatter plot +close all + +f1 = figure('Position',[200 500 1200 800],'DefaultAxesFontSize',12,'renderer','painters'); +t = tiledlayout(2,3,'TileSpacing','tight','Padding','tight'); + +s1=nexttile(1); +hold on +l1=plot(socBoth(:,1),yLs,'-r','LineWidth',2); +l2=plot(socBoth(:,1),yLs+dLs,'--r','LineWidth',1); +plot(socBoth(:,1),yLs-dLs,'--r','LineWidth',1); +l3=plot(socBoth(:,1),yCs,'-g','LineWidth',2); +l4=plot(socBoth(:,1),yCs+dCs,'--g','LineWidth',1); +plot(socBoth(:,1),yCs-dCs,'--g','LineWidth',1); +scatter(socBoth(:,1),socBoth(:,2),60,'filled','MarkerFaceColor','k','MarkerEdgeColor','w'); +xlim([min(socBoth(:,1))-1,max(socBoth(:,1))+1]); +ylim([min(socBoth(:,2))-1,max(socBoth(:,2))+6]); + +legend([l1,l2,l3,l4],{['y=',num2str(pLs(1)),'x+',num2str(pLs(2))], ... + ['Mean standard error ',num2str(mean(dLs))], ... + ['y=',num2str(pCs(1)),'x^2+',num2str(pCs(2)),'x+',num2str(pCs(3))], ... + ['Mean standard error ',num2str(mean(dCs))]}); + +xlabel('Aircraft speed (m s^{-1})') +ylabel('Peak at number of non-zeros') + +title(['SOCRATES binned']); + +grid on +box on + +s2=nexttile(2); +hold on +l1=plot(otBoth(:,1),yLo,'-r','LineWidth',2); +l2=plot(otBoth(:,1),yLo+dLo,'--r','LineWidth',1); +plot(otBoth(:,1),yLo-dLo,'--r','LineWidth',1); +l3=plot(otBoth(:,1),yCo,'-g','LineWidth',2); +l4=plot(otBoth(:,1),yCo+dCo,'--g','LineWidth',1); +plot(otBoth(:,1),yCo-dCo,'--g','LineWidth',1); +scatter(otBoth(:,1),otBoth(:,2),60,'filled','MarkerFaceColor','k','MarkerEdgeColor','w'); +xlim([min(otBoth(:,1))-1,max(otBoth(:,1))+1]); +ylim([min(otBoth(:,2))-1,max(otBoth(:,2))+6]); + +legend([l1,l2,l3,l4],{['y=',num2str(pLo(1)),'x+',num2str(pLo(2))], ... + ['Mean standard error ',num2str(mean(dLo))], ... + ['y=',num2str(pCo(1)),'x^2+',num2str(pCo(2)),'x+',num2str(pCo(3))], ... + ['Mean standard error ',num2str(mean(dCo))]}); + +xlabel('Aircraft speed (m s^{-1})') +ylabel('Peak at number of non-zeros') + +title(['OTREC binned']); + +grid on +box on + +s3=nexttile(3); +hold on +l1=plot(comb(:,1),yLb,'-r','LineWidth',2); +l2=plot(comb(:,1),yLb+dLb,'--r','LineWidth',1); +plot(comb(:,1),yLb-dLb,'--r','LineWidth',1); +l3=plot(comb(:,1),yCb,'-g','LineWidth',2); +l4=plot(comb(:,1),yCb+dCb,'--g','LineWidth',1); +plot(comb(:,1),yCb-dCb,'--g','LineWidth',1); +scatter(comb(:,1),comb(:,2),60,'filled','MarkerFaceColor','k','MarkerEdgeColor','w'); +xlim([min(comb(:,1))-1,max(comb(:,1))+1]); +ylim([min(comb(:,2))-1,max(comb(:,2))+6]); + +legend([l1,l2,l3,l4],{['y=',num2str(pLb(1)),'x+',num2str(pLb(2))], ... + ['Mean standard error ',num2str(mean(dLb))], ... + ['y=',num2str(pCb(1)),'x^2+',num2str(pCb(2)),'x+',num2str(pCb(3))], ... + ['Mean standard error ',num2str(mean(dCb))]}); + +xlabel('Aircraft speed (m s^{-1})') +ylabel('Peak at number of non-zeros') + +title(['SOCRATES and OTREC binned']); + +grid on +box on + +%% Cases +socrates=load('/scr/virga2/rsfdata/projects/socrates/hcr/qc3/cfradial/v3.2_full/spectralMoments/smoothingTest/socrates_smoothingAnalysis_everyOther_aircraftSpeed_cases.mat'); +otrec=load('/scr/sleet2/rsfdata/projects/otrec/hcr/qc3/cfradial/v3.2_full/spectralMoments/smoothingTest/otrec_smoothingAnalysis_everyOther_aircraftSpeed_cases.mat'); + +%% Combine +socBoth=cat(2,socrates.velAircMeanC,socrates.peakNumCases); +socBoth(any(isnan(socBoth),2),:)=[]; +otBoth=cat(2,otrec.velAircMeanC,otrec.peakNumCases); +otBoth(any(isnan(otBoth),2),:)=[]; + +comb=cat(1,socBoth,otBoth); +comb=sortrows(comb); + + +%% Fits + +[pLs,sLs]=polyfit(socBoth(:,1),socBoth(:,2),1); +[yLs,dLs]=polyval(pLs,socBoth(:,1),sLs); + +[pCs,sCs]=polyfit(socBoth(:,1),socBoth(:,2),2); +[yCs,dCs]=polyval(pCs,socBoth(:,1),sCs); + +[pLo,sLo]=polyfit(otBoth(:,1),otBoth(:,2),1); +[yLo,dLo]=polyval(pLo,otBoth(:,1),sLo); + +[pCo,sCo]=polyfit(otBoth(:,1),otBoth(:,2),2); +[yCo,dCo]=polyval(pCo,otBoth(:,1),sCo); + +[pLb,sLb]=polyfit(comb(:,1),comb(:,2),1); +[yLb,dLb]=polyval(pLb,comb(:,1),sLb); + +[pCb,sCb]=polyfit(comb(:,1),comb(:,2),2); +[yCb,dCb]=polyval(pCb,comb(:,1),sCb); + +%% Scatter plot + +s4=nexttile(4); +hold on +l1=plot(socBoth(:,1),yLs,'-r','LineWidth',2); +l2=plot(socBoth(:,1),yLs+dLs,'--r','LineWidth',1); +plot(socBoth(:,1),yLs-dLs,'--r','LineWidth',1); +l3=plot(socBoth(:,1),yCs,'-g','LineWidth',2); +l4=plot(socBoth(:,1),yCs+dCs,'--g','LineWidth',1); +plot(socBoth(:,1),yCs-dCs,'--g','LineWidth',1); +scatter(socBoth(:,1),socBoth(:,2),60,'filled','MarkerFaceColor','k','MarkerEdgeColor','w'); +xlim([min(socBoth(:,1))-1,max(socBoth(:,1))+1]); +ylim([min(socBoth(:,2))-1,max(socBoth(:,2))+6]); + +legend([l1,l2,l3,l4],{['y=',num2str(pLs(1)),'x+',num2str(pLs(2))], ... + ['Mean standard error ',num2str(mean(dLs))], ... + ['y=',num2str(pCs(1)),'x^2+',num2str(pCs(2)),'x+',num2str(pCs(3))], ... + ['Mean standard error ',num2str(mean(dCs))]}); + +xlabel('Aircraft speed (m s^{-1})') +ylabel('Peak at number of non-zeros') + +title(['SOCRATES cases']); + +grid on +box on + +s5=nexttile(5); +hold on +l1=plot(otBoth(:,1),yLo,'-r','LineWidth',2); +l2=plot(otBoth(:,1),yLo+dLo,'--r','LineWidth',1); +plot(otBoth(:,1),yLo-dLo,'--r','LineWidth',1); +l3=plot(otBoth(:,1),yCo,'-g','LineWidth',2); +l4=plot(otBoth(:,1),yCo+dCo,'--g','LineWidth',1); +plot(otBoth(:,1),yCo-dCo,'--g','LineWidth',1); +scatter(otBoth(:,1),otBoth(:,2),60,'filled','MarkerFaceColor','k','MarkerEdgeColor','w'); +xlim([min(otBoth(:,1))-1,max(otBoth(:,1))+1]); +ylim([min(otBoth(:,2))-1,max(otBoth(:,2))+6]); + +legend([l1,l2,l3,l4],{['y=',num2str(pLo(1)),'x+',num2str(pLo(2))], ... + ['Mean standard error ',num2str(mean(dLo))], ... + ['y=',num2str(pCo(1)),'x^2+',num2str(pCo(2)),'x+',num2str(pCo(3))], ... + ['Mean standard error ',num2str(mean(dCo))]}); + +xlabel('Aircraft speed (m s^{-1})') +ylabel('Peak at number of non-zeros') + +title(['OTREC cases']); + +grid on +box on + +s6=nexttile(6); +hold on +l1=plot(comb(:,1),yLb,'-r','LineWidth',2); +l2=plot(comb(:,1),yLb+dLb,'--r','LineWidth',1); +plot(comb(:,1),yLb-dLb,'--r','LineWidth',1); +l3=plot(comb(:,1),yCb,'-g','LineWidth',2); +l4=plot(comb(:,1),yCb+dCb,'--g','LineWidth',1); +plot(comb(:,1),yCb-dCb,'--g','LineWidth',1); +scatter(comb(:,1),comb(:,2),60,'filled','MarkerFaceColor','k','MarkerEdgeColor','w'); +xlim([min(comb(:,1))-1,max(comb(:,1))+1]); +ylim([min(comb(:,2))-1,max(comb(:,2))+6]); + +legend([l1,l2,l3,l4],{['y=',num2str(pLb(1)),'x+',num2str(pLb(2))], ... + ['Mean standard error ',num2str(mean(dLb))], ... + ['y=',num2str(pCb(1)),'x^2+',num2str(pCb(2)),'x+',num2str(pCb(3))], ... + ['Mean standard error ',num2str(mean(dCb))]}); + +xlabel('Aircraft speed (m s^{-1})') +ylabel('Peak at number of non-zeros') + +title(['SOCRATES and OTREC cases']); + +grid on +box on + +set(gcf,'PaperPositionMode','auto') +print(f1,[figdir,'smoothingAnalysis_everyOther_aircraftSpeed'],'-dpng','-r0');