diff --git a/Misc/NormalizationCheck/LIGOnoise.m b/Misc/NormalizationCheck/LIGOnoise.m new file mode 100644 index 0000000..aae5ba7 --- /dev/null +++ b/Misc/NormalizationCheck/LIGOnoise.m @@ -0,0 +1,63 @@ +function [outNoise, PSD] = LIGOnoise(N, Fs, noise_num, noisefile) +%Function to create colored noise using LIGO Design Sensitivities +% Design PSD is modified between 15 Hz and 700Hz. +% Input: N = Total number of samples, +% Fs = Sampling Frequency, +% (Optional) noise_num = noise realization number from a pre-created noise realizations file +% (Optional) noisefile = pre-created noise realizations filename +% Output: outNoise = colored noise vector, +% PSD = two-sided PSD vector for positive DFT frequencies + +% Raghav Girgaonkar, April 2023 + +%Load PSD +y = load('iLIGOSensitivity.txt','-ascii'); +% freqs = y(:,1); +% sqrtPSD = y(:,2); + +%Turn one-sided sensitivity to two-sided +y(:,2) = (1/sqrt(2))*y(:,2); + +% Interpolate sensitivity curve to positive DFT frequencies +minF = min(y(:,1)); +maxF = max(y(:,1)); +if minF ~= 0 +% f=0 does not exist, put it in +y = [0, y(1,2);... + y]; +end +if maxF < Fs/2 + error('High frequency limit requested is higher than supplied'); +end + + +%Positive DFT frequencies +kNyq = floor(N/2)+1; +fvec = (0:(kNyq-1))*Fs/N; + +%% Interpolation +interPSD = interp1(y(:,1),y(:,2), fvec); + +%% Modifications, change cutoff frequencies as needed +minidx = find(fvec<=30, 1, 'last' ); +maxidx = find(fvec<=700, 1, 'last' ); + +Sn50 = interPSD(minidx); +Sn700 = interPSD(maxidx); + +interPSD(1:minidx) = Sn50; +interPSD(maxidx:end) = Sn700; + +%Convert Amplitude spectral density to Power Spectral Density +PSD = interPSD.^2; + +%% Make colored Noise +fltrOrdr = 10000; + +outNoise_t = statgaussnoisegen(N,[fvec(:),PSD(:)],fltrOrdr,Fs, noise_num, noisefile); + +outNoise = outNoise_t(fltrOrdr+1:end - fltrOrdr); + + + + diff --git a/Misc/NormalizationCheck/allparamfiles.json b/Misc/NormalizationCheck/allparamfiles.json new file mode 100644 index 0000000..a493762 --- /dev/null +++ b/Misc/NormalizationCheck/allparamfiles.json @@ -0,0 +1,5 @@ +{ +"signalparamfile": "signal.json", +"psoparamfile": "pso.json", +"filenames": "files.json" +} diff --git a/Misc/NormalizationCheck/animate_matchedfiltering.m b/Misc/NormalizationCheck/animate_matchedfiltering.m new file mode 100644 index 0000000..945e026 --- /dev/null +++ b/Misc/NormalizationCheck/animate_matchedfiltering.m @@ -0,0 +1,127 @@ +%% Script to create an animation of matched filtering timeseries. + +%% Define Data Parameters +sampFreq = 4096; +datalen = 128; +N = datalen*sampFreq; +frange = [30,700]; +fpos = (0:floor(N/2))*(1/datalen); +timeVec = (0:(datalen*sampFreq - 1))/sampFreq; + +%% Create Colored Gaussian Noise +[noise, PSD] = LIGOnoise(N, sampFreq, 1, 'sample'); + + +%Create Total PSD and Transfer Function +negFStrt = 1-mod(N,2); +kNyq = floor(N/2)+1; + +PSDtotal = [PSD, PSD((kNyq-negFStrt):-1:2)]; +TFtotal = 1./sqrt(PSDtotal); + + +%Pre-calculate frequency magnitude and phase difference terms + [A,avec, phaseDiff] = preprocessing(frange(1),frange(2),fpos, datalen, N); + +%Create general normalization factor N_fW +AbysqrtPSD = A.*TFtotal; +innProd = (1/N)*(AbysqrtPSD)*AbysqrtPSD'; +genNormfacSqr = real(innProd); +N_fW = 1/sqrt(genNormfacSqr); + +%% Generate Signal to Inject + chirptimes = [29.6373, 1.1045]; + ta = 34; + initial_phase = 0; + phase = 0; + snr = 12; + +signal = gen2PNtemplate_tau(fpos, ta, phase, frange(1), frange(2),chirptimes(1), chirptimes(2),datalen,initial_phase,snr, N, A, avec, PSDtotal); + +signal = signal*sqrt(sampFreq); + +figure; +plot(timeVec, noise); +hold on; +plot(timeVec, 30*signal); +xlabel('Time (s)'); +ylabel('Strain h(t)'); +axis tight; +ax = gca(); ax.XAxis.FontSize = 20; ax.YAxis.FontSize = 20; +%Save whitened and normalized signal and noise +%Window first +whtndsignal = ifft((1/sqrt(sampFreq))*fft(signal).*TFtotal); + +%Inject signal into data +strain = noise + signal; + +%Whiten and normalize strain +rolloff = 0.5; +fmin = frange(1); +strainwin = strain.*tukeywin(length(strain),rolloff*sampFreq/length(strain))'; +whtndstrain = ifft((1/sqrt(sampFreq))*fft(strainwin).*TFtotal); + + +%Make GIF +% Define the filename for the gif +giffilename = 'matchedfiltering_animation.gif'; +time_lapse = 0.0001; +title_text = 'Generation of Matched-Filtering Timeseries'; + +% Create a new figure window +fig = figure; +mftimeseries = zeros(1,N); +for i = 1:1*sampFreq:N + q0 = gen2PNtemplate_tau(fpos, i/sampFreq, 0, frange(1), frange(2),chirptimes(1), chirptimes(2),datalen,initial_phase,1, N, A, avec, PSDtotal); + whtndq0 = ifft(fft(q0).*TFtotal); + q1 = gen2PNtemplate_tau(fpos, i/sampFreq, pi/2, frange(1), frange(2),chirptimes(1), chirptimes(2),datalen,initial_phase,1, N, A, avec, PSDtotal); + whtndq1 = ifft(fft(q1).*TFtotal); + +% I0 = innerproduct(strain, q0, PSDtotal, sampFreq); +% I1 = innerproduct(strain, q1, PSDtotal, sampFreq); + + I0 = real((1/N)*sum(fft(whtndstrain).*conj(fft(whtndq0)))); + I1 = real((1/N)*sum(fft(whtndstrain).*conj(fft(whtndq1)))); + + I = sqrt(I0^2 + I1^2); + + mftimeseries(i) = I; + + + hold on; + plot(timeVec, whtndstrain); hold on; + plot(timeVec, 5*whtndsignal, Color='red' ); hold on; + plot(timeVec, 20*whtndq0,Color='yellow'); hold on; + plot(timeVec(1:i), mftimeseries(1:i), Color='black',LineWidth=1.5); hold on; + ylim([-4,12]); + xlim([0,128]); + xlabel("Time (s)"); + ylabel("Whitened Strain h(t)"); + ax = gca(); ax.XAxis.FontSize = 20; ax.YAxis.FontSize = 20; +% axis tight; + + title(title_text,'FontSize',16,'Color',[0 0 0]), hold on; + + drawnow; + %Capture frame + % Capture the current frame and save it to the gif + frame = getframe(fig); + im = frame2im(frame); + [imind,cm] = rgb2ind(im,256); + + if i == 1 + imwrite(imind,cm,giffilename,'gif', 'Loopcount',Inf,'DelayTime',time_lapse); + else + imwrite(imind,cm,giffilename,'gif','WriteMode','append','DelayTime',time_lapse); + end + + clf(fig); + +end + +% +% figure; +% plot(timeVec, noise); +% hold on; +% plot(timeVec, 10*signal); +% axis tight; diff --git a/Misc/NormalizationCheck/animate_matchedfiltering_glitch.m b/Misc/NormalizationCheck/animate_matchedfiltering_glitch.m new file mode 100644 index 0000000..40edb56 --- /dev/null +++ b/Misc/NormalizationCheck/animate_matchedfiltering_glitch.m @@ -0,0 +1,128 @@ +%% Script to create an animation of matched filtering timeseries. + +%% Define Data Parameters +sampFreq = 4096; +datalen = 128; +N = datalen*sampFreq; +frange = [30,700]; +fpos = (0:floor(N/2))*(1/datalen); +timeVec = (0:(datalen*sampFreq - 1))/sampFreq; + +%% Create Colored Gaussian Noise +[noise, PSD] = LIGOnoise(N, sampFreq, 1, 'sample'); + + +%Create Total PSD and Transfer Function +negFStrt = 1-mod(N,2); +kNyq = floor(N/2)+1; + +PSDtotal = [PSD, PSD((kNyq-negFStrt):-1:2)]; +TFtotal = 1./sqrt(PSDtotal); + + +%Pre-calculate frequency magnitude and phase difference terms + [A,avec, phaseDiff] = preprocessing(frange(1),frange(2),fpos, datalen, N); + +%Create general normalization factor N_fW +AbysqrtPSD = A.*TFtotal; +innProd = (1/N)*(AbysqrtPSD)*AbysqrtPSD'; +genNormfacSqr = real(innProd); +N_fW = 1/sqrt(genNormfacSqr); + +%% Generate Signal to Inject + chirptimes = [0.001,0.01]; + qchirptimes = [11.6373, 1.1045]; + ta = 54; + initial_phase = 0; + phase = 0; + snr = 360; + +signal = gen2PNtemplate_tau_negative(fpos, ta, phase, frange(1), frange(2),chirptimes(1), chirptimes(2),datalen,initial_phase,snr, N, A, avec, PSDtotal); + +signal = signal*sqrt(sampFreq); + +% figure; +% plot(timeVec, noise); +% hold on; +% plot(timeVec, 30*signal); +% xlabel('Time (s)'); +% ylabel('Strain h(t)'); +% axis tight; +% ax = gca(); ax.XAxis.FontSize = 20; ax.YAxis.FontSize = 20; +%Save whitened and normalized signal and noise +%Window first +whtndsignal = ifft((1/sqrt(sampFreq))*fft(signal).*TFtotal); + +%Inject signal into data +strain = noise + signal; + +%Whiten and normalize strain +rolloff = 0.5; +fmin = frange(1); +strainwin = strain.*tukeywin(length(strain),rolloff*sampFreq/length(strain))'; +whtndstrain = ifft((1/sqrt(sampFreq))*fft(strainwin).*TFtotal); + + +%Make GIF +% Define the filename for the gif +giffilename = 'matchedfiltering_glitch_animation.gif'; +time_lapse = 0.0001; +title_text = 'Generation of Matched-Filtering Timeseries'; + +% Create a new figure window +fig = figure; +mftimeseries = zeros(1,N); +for i = 1:1*sampFreq:N + q0 = gen2PNtemplate_tau(fpos, i/sampFreq, 0, frange(1), frange(2),qchirptimes(1), qchirptimes(2),datalen,initial_phase,1, N, A, avec, PSDtotal); + whtndq0 = ifft(fft(q0).*TFtotal); + q1 = gen2PNtemplate_tau(fpos, i/sampFreq, pi/2, frange(1), frange(2),qchirptimes(1), qchirptimes(2),datalen,initial_phase,1, N, A, avec, PSDtotal); + whtndq1 = ifft(fft(q1).*TFtotal); + +% I0 = innerproduct(strain, q0, PSDtotal, sampFreq); +% I1 = innerproduct(strain, q1, PSDtotal, sampFreq); + + I0 = real((1/N)*sum(fft(whtndstrain).*conj(fft(whtndq0)))); + I1 = real((1/N)*sum(fft(whtndstrain).*conj(fft(whtndq1)))); + + I = sqrt(I0^2 + I1^2); + + mftimeseries(i) = I; + + + hold on; + plot(timeVec, whtndstrain); hold on; + plot(timeVec, whtndsignal, Color='red' ); hold on; + plot(timeVec, 20*whtndq0,Color='yellow'); hold on; + plot(timeVec(1:i), mftimeseries(1:i), Color='black',LineWidth=1.5); hold on; + ylim([-4,15]); + xlim([0,128]); + xlabel("Time (s)"); + ylabel("Whitened Strain h(t)"); + ax = gca(); ax.XAxis.FontSize = 20; ax.YAxis.FontSize = 20; +% axis tight; + + title(title_text,'FontSize',16,'Color',[0 0 0]), hold on; + + drawnow; + %Capture frame + % Capture the current frame and save it to the gif + frame = getframe(fig); + im = frame2im(frame); + [imind,cm] = rgb2ind(im,256); + + if i == 1 + imwrite(imind,cm,giffilename,'gif', 'Loopcount',Inf,'DelayTime',time_lapse); + else + imwrite(imind,cm,giffilename,'gif','WriteMode','append','DelayTime',time_lapse); + end + + clf(fig); + +end + +% +% figure; +% plot(timeVec, noise); +% hold on; +% plot(timeVec, 10*signal); +% axis tight; diff --git a/Misc/NormalizationCheck/boundary_plot.m b/Misc/NormalizationCheck/boundary_plot.m new file mode 100644 index 0000000..b268220 --- /dev/null +++ b/Misc/NormalizationCheck/boundary_plot.m @@ -0,0 +1,78 @@ +function [boundary_fig] = boundary_plot +% Creates and returns a Boundary Plot as a figure in the Chirp-time (Tau1.5 and Tau0) space +% Search range of component masses is 1.4 to 30 Solar Masses + +%Constants +c = 3*10^8; +Msolar = 1.989*10^30; +G = 6.6743*10^-11; +%Frequency min max values +fmin = 30; +fmax = 700; +%Mass vectors +m_lin = linspace(1.4,30,10000); +m1 = m_lin; +m2 = m1; +%Calculate Mass Terms +m1 = m1*Msolar; +m2 = m2*Msolar; +M = (m1 + m2); +u = m1.*m2./M; +n = u./M; +%Calculate Chirp Times firt when m1 = m2 +tau0_1 = (5/(256*pi))*(1/fmin)*((G*M*pi*fmin/c^3).^(-5/3)).*(1./n); + +% tau1 = (5/(192*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-1))*(1/n)*((743/336)+ (11*n/4)); + +tau1p5_1 = (1/8)*(1/fmin)*((G*M*pi*fmin/c^3).^(-2/3)).*(1./n); + +% tau2 = (5/(128*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-1/3))*(1/n)*((3058673/1016064) + (5429*n/1008) + (617*n*n/144)); + +%Boundary Array for m2 = 1.4 +%Mass vectors +m1 = m_lin; +m2 = m1(1); +%Calculate Mass Terms +m1 = m1*Msolar; +m2 = m2*Msolar; +M = (m1 + m2); +u = (m1*m2)./M; +n = u./M; +tau0_2 = (5/(256*pi))*(1/fmin)*((G*M*pi*fmin/c^3).^(-5/3)).*(1./n); +tau1p5_2 = (1/8)*(1/fmin)*((G*M*pi*fmin/c^3).^(-2/3)).*(1./n); + +%Boundary Array for m1 = 30 +%Mass vectors +m2 = m_lin; +m1 = m2(end); +%Calculate Mass Terms +m1 = m1*Msolar; +m2 = m2*Msolar; +M = (m1 + m2); +u = (m1*m2)./M; +n = u./M; +tau0_3 = (5/(256*pi))*(1/fmin)*((G*M*pi*fmin/c^3).^(-5/3)).*(1./n); +tau1p5_3 = (1/8)*(1/fmin)*((G*M*pi*fmin/c^3).^(-2/3)).*(1./n); + +%Final Tau0 and Tau1.5 arrays + +Tau0 = [tau0_1, tau0_2, tau0_3]; + +Tau1p5 = [tau1p5_1, tau1p5_2, tau1p5_3]; + +%Scatterplot +sz = 5; +c = 'black'; +% boundary_fig = figure; +% hold on; +% t1 = scatter(tau0_1, tau1p5_1,sz, 'red', 'filled', 'HandleVisibility','off'); +% hold on; +% t2 = scatter(tau0_2, tau1p5_2,sz, c, 'filled', 'HandleVisibility','off'); +% hold on; +% t3 = scatter(tau0_3, tau1p5_3,sz, c, 'filled', 'HandleVisibility','off'); +% hold off; +boundary_fig = scatter(Tau0, Tau1p5,sz,c,"filled", 'DisplayName','Tau Space Boundary'); +% xlabel('\tau_0'); +% ylabel('\tau_{1.5}'); +% title("Boundary of Region of Physical Parameters in \tau_0 and \tau_{1.5} Space"); +% saveas(gcf,"boundary-tauspace.pdf"); \ No newline at end of file diff --git a/Misc/NormalizationCheck/crcbgwpso_mass.m b/Misc/NormalizationCheck/crcbgwpso_mass.m new file mode 100644 index 0000000..878d6ad --- /dev/null +++ b/Misc/NormalizationCheck/crcbgwpso_mass.m @@ -0,0 +1,121 @@ +function outResults = crcbgwpso_mass(inParams,psoParams,nRuns, sampling_freq) +%Regression of 2PNWaveform using Mass Space PSO +% Input: inParams: Struct containing data and signal parameters +% psoParams: Struct containing PSO parameters +% nRuns: Number of PSO iterations +% sampling_freq: Sampling Frequency of data +% Output: outResults: Struct containing the following parameters +% 'allRunsOutput': An N element struct array containing results from each PSO +% run. The fields of this struct are: +% 'fitVal': The fitness value. +% 'qcCoefs': The coefficients [m1, m2]. +% 'estSig': The estimated signal. +% 'totalFuncEvals': The total number of fitness +% evaluations. +% 'bestRun': The best run. +% 'bestFitness': best fitness from the best run. +% 'bestSig' : The signal estimated in the best run. +% 'bestQcCoefs' : [m1, m2] found in the best run. + +%Raghav Girgaonkar, April 2023 + + +nSamples = length(inParams.dataX); + +fHandle = @(x) psofitfunc(x,inParams); + +params = inParams; + +nDim = 2; +outStruct = struct('bestLocation',[],... + 'bestFitness', [],... + 'totalFuncEvals',[],... + 'allBestFit',[],... + 'allBestLoc',[]); + +outResults = struct('allRunsOutput',struct('fitVal', [],... + 'qcCoefs',zeros(1,2),... + 'estTa',[],... + 'estSig',zeros(1,nSamples),... + 'totalFuncEvals',[],... + 'allBestFit',zeros(1,psoParams.maxSteps),... + 'allBestLoc',zeros(nDim,psoParams.maxSteps)),... + 'bestRun',[],... + 'bestFitness',[],... + 'bestSig', zeros(1,nSamples),... + 'bestQcCoefs',zeros(1,2),... + 'estAmp',[],... + 'estPhase',[],... + 'bestAmp',[],... + 'bestPhase',[],... + 'bestTime',[]); + +%Allocate storage for outputs: results from all runs are stored +for lpruns = 1:nRuns + outStruct(lpruns) = outStruct(1); + outResults.allRunsOutput(lpruns) = outResults.allRunsOutput(1); +end +%Independent runs of PSO in parallel. Change 'parfor' to 'for' if the +%parallel computing toolbox is not available. +% fprintf("Running PSO\n"); +% parpool(nruns); +parfor lpruns = 1:nRuns + %Reset random number generator for each worker + rng(lpruns); + outStruct(lpruns)=crcbpso(fHandle,nDim,psoParams,2); +end + +%Prepare output +fitVal = zeros(1,nRuns); +for lpruns = 1:nRuns + outResults.allRunsOutput(lpruns).allBestFit = outStruct(lpruns).allBestFit; + outResults.allRunsOutput(lpruns).allBestLoc = outStruct(lpruns).allBestLoc; + fitVal(lpruns) = outStruct(lpruns).bestFitness; + outResults.allRunsOutput(lpruns).fitVal = fitVal(lpruns); + [~,masses,~] = fHandle(outStruct(lpruns).bestLocation); + + outResults.allRunsOutput(lpruns).qcCoefs = masses; + + %Get parameter estimates + [estAmp, estTa, estPhase] = getparamestimates_mass(chirptimes, params); + %Calculate time using sampling freq and ta_index +% estTa = ta_index/sampling_freq; + + outResults.allRunsOutput(lpruns).estTa = estTa; + m1 = masses(1); + m2 = masses(2); +% phaseq0 = gen2PNwaveform(params.fpos, estTa, 0, params.frange(1), params.frange(2), m1,... +% m2,params.datalen,0,1,params.N,params.avec, params.normfac); +% +% fftq0 = phaseq0; +% fftq1 = phaseq0.*params.phaseDiff; +% +% +% yq0 = innerprodpsd(fftq0, params.fftdataYbyPSD); +% yq1 = innerprodpsd(fftq1, params.fftdataYbyPSD); +% estPhase = atan2(yq1,yq0); + outResults.allRunsOutput(lpruns).estPhase = estPhase; + + %Estimated Amplitude +% estAmp = cos(estPhase)*yq0 + sin(estPhase)*yq1; + outResults.allRunsOutput(lpruns).estAmp = estAmp; + %Estimated Signal + + estSigphase = gen2PNwaveform(params.fpos, estTa, estPhase, params.frange(1), params.frange(2), m1,... + m2,params.datalen,0,estAmp,params.N,params.avec, params.normfac); + estSigfourier = (params.A).*estSigphase; + estSig = ifft(estSigfourier); + + outResults.allRunsOutput(lpruns).estSig = estSig; + outResults.allRunsOutput(lpruns).totalFuncEvals = outStruct(lpruns).totalFuncEvals; +end +%Find the best run +[~,bestRun] = min(fitVal(:)); +outResults.bestRun = bestRun; +outResults.bestFitness = outResults.allRunsOutput(bestRun).fitVal; +outResults.bestSig = outResults.allRunsOutput(bestRun).estSig; +outResults.bestAmp = outResults.allRunsOutput(bestRun).estAmp; +outResults.bestPhase = outResults.allRunsOutput(bestRun).estPhase; +outResults.bestQcCoefs = outResults.allRunsOutput(bestRun).qcCoefs; +outResults.bestTime = outResults.allRunsOutput(bestRun).estTa; + diff --git a/Misc/NormalizationCheck/crcbgwpso_tau.m b/Misc/NormalizationCheck/crcbgwpso_tau.m new file mode 100644 index 0000000..02238dd --- /dev/null +++ b/Misc/NormalizationCheck/crcbgwpso_tau.m @@ -0,0 +1,124 @@ +function outResults = crcbgwpso_tau(inParams,psoParams,nRuns, sampling_freq) +%Regression of 2PNWaveform using Chirp-Time Space PSO +% Input: inParams: Struct containing data and signal parameters +% psoParams: Struct containing PSO parameters +% nRuns: Number of PSO iterations +% sampling_freq: Sampling Frequency of data +% Output: outResults: Struct containing the following parameters +% 'allRunsOutput': An N element struct array containing results from each PSO +% run. The fields of this struct are: +% 'fitVal': The fitness value. +% 'qcCoefs': The coefficients [tau0, tau1.5]. +% 'estSig': The estimated signal. +% 'totalFuncEvals': The total number of fitness +% evaluations. +% 'bestRun': The best run. +% 'bestFitness': best fitness from the best run. +% 'bestSig' : The signal estimated in the best run. +% 'bestQcCoefs' : [tau0, tau1.5] found in the best run. + +%Raghav Girgaonkar, April 2023 + +nSamples = length(inParams.dataX); + +fHandle = @(x) psofitfunc_tau(x,inParams); + +params = inParams; + +nDim = 2; +outStruct = struct('bestLocation',[],... + 'bestFitness', [],... + 'totalFuncEvals',[],... + 'allBestFit',[],... + 'allBestLoc',[]); + +outResults = struct('allRunsOutput',struct('fitVal', [],... + 'qcCoefs',zeros(1,2),... + 'estTa',[],... + 'estSig',zeros(1,nSamples),... + 'totalFuncEvals',[],... + 'allBestFit',zeros(1,psoParams.maxSteps),... + 'allBestLoc',zeros(nDim,psoParams.maxSteps)),... + 'bestRun',[],... + 'bestFitness',[],... + 'bestSig', zeros(1,nSamples),... + 'bestQcCoefs',zeros(1,2),... + 'estAmp',[],... + 'estPhase',[],... + 'bestAmp',[],... + 'bestPhase',[],... + 'bestTime',[]); + +%Allocate storage for outputs: results from all runs are stored +for lpruns = 1:nRuns + outStruct(lpruns) = outStruct(1); + outResults.allRunsOutput(lpruns) = outResults.allRunsOutput(1); +end +%Independent runs of PSO in parallel. Change 'parfor' to 'for' if the +%parallel computing toolbox is not available. +% fprintf("Running PSO\n"); +parpool(nRuns); +parfor lpruns = 1:nRuns + %Reset random number generator for each worker + rng(lpruns); + outStruct(lpruns)=crcbpso(fHandle,nDim,psoParams,2); +end + +%Prepare output +fitVal = zeros(1,nRuns); +% sampling_interval = 1/sampling_freq; +%Time vectors for signal and total series +% timeVecSig = (0:(sampling_freq*T_sig - 1))/sampling_freq; +for lpruns = 1:nRuns + outResults.allRunsOutput(lpruns).allBestFit = outStruct(lpruns).allBestFit; + outResults.allRunsOutput(lpruns).allBestLoc = outStruct(lpruns).allBestLoc; + fitVal(lpruns) = outStruct(lpruns).bestFitness; + outResults.allRunsOutput(lpruns).fitVal = fitVal(lpruns); + [~,chirptimes,~] = fHandle(outStruct(lpruns).bestLocation); +% Q = qcCoefs +% index = ta_index + [estAmp, estTa, estPhase] = getparamestimates(chirptimes, params); + outResults.allRunsOutput(lpruns).qcCoefs = chirptimes; + %Calculate time using sampling freq and ta_index +% estTa = ta_index/sampling_freq; + + outResults.allRunsOutput(lpruns).estTa = estTa; + tau0 = chirptimes(1); + tau1p5 = chirptimes(2); +% phaseq0 = gen2PNwaveform_tau(params.fpos, estTa, 0, params.frange(1), params.frange(2), tau0,... +% tau1p5,params.datalen,0,1,params.N,params.avec, params.normfac); +% fftq0 = phaseq0; +% fftq1 = phaseq0.*params.phaseDiff; +% %Estimated Phase +% % yq0 = inParams.dataY*q0(:); +% % yq1 = inParams.dataY*q1(:); +% yq0 = innerprodpsd(fftq0, params.fftdataYbyPSD); +% yq1 = innerprodpsd(fftq1, params.fftdataYbyPSD); +% estPhase = atan2(yq1,yq0); + outResults.allRunsOutput(lpruns).estPhase = estPhase; + + %Estimated Amplitude +% estAmp = cos(estPhase)*yq0 + sin(estPhase)*yq1; + outResults.allRunsOutput(lpruns).estAmp = estAmp; + + %Estimated Signal +% estSigTemp = genqc(timeVecSig,1,qcCoefs,estPhase); + estSigphase = gen2PNtemplate_tau(params.fpos, estTa, estPhase, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,estAmp,params.N,params.A,params.avec, params.PSDtotal); + estSigfourier = (params.A).*estSigphase; + estSig = ifft(estSigfourier); +% estSigTemp_shifted = [zeros(1,floor(estTa*sampling_freq)-1), estSigTemp, zeros(1, nSamples - nSamplesSig - floor(estTa*sampling_freq)+1)]; +% estSig = estAmp*estSigTemp; + outResults.allRunsOutput(lpruns).estSig = estSig; + outResults.allRunsOutput(lpruns).totalFuncEvals = outStruct(lpruns).totalFuncEvals; +end +%Find the best run +[~,bestRun] = min(fitVal(:)); +outResults.bestRun = bestRun; +outResults.bestFitness = outResults.allRunsOutput(bestRun).fitVal; +outResults.bestSig = outResults.allRunsOutput(bestRun).estSig; +outResults.bestAmp = outResults.allRunsOutput(bestRun).estAmp; +outResults.bestPhase = outResults.allRunsOutput(bestRun).estPhase; +outResults.bestQcCoefs = outResults.allRunsOutput(bestRun).qcCoefs; +outResults.bestTime = outResults.allRunsOutput(bestRun).estTa; + diff --git a/Misc/NormalizationCheck/crcbgwpso_tau_negative.m b/Misc/NormalizationCheck/crcbgwpso_tau_negative.m new file mode 100644 index 0000000..d5c47fa --- /dev/null +++ b/Misc/NormalizationCheck/crcbgwpso_tau_negative.m @@ -0,0 +1,124 @@ +function outResults = crcbgwpso_tau_negative(inParams,psoParams,nRuns, sampling_freq) +%Regression of 2PNWaveform using Chirp-Time Space PSO +% Input: inParams: Struct containing data and signal parameters +% psoParams: Struct containing PSO parameters +% nRuns: Number of PSO iterations +% sampling_freq: Sampling Frequency of data +% Output: outResults: Struct containing the following parameters +% 'allRunsOutput': An N element struct array containing results from each PSO +% run. The fields of this struct are: +% 'fitVal': The fitness value. +% 'qcCoefs': The coefficients [tau0, tau1.5]. +% 'estSig': The estimated signal. +% 'totalFuncEvals': The total number of fitness +% evaluations. +% 'bestRun': The best run. +% 'bestFitness': best fitness from the best run. +% 'bestSig' : The signal estimated in the best run. +% 'bestQcCoefs' : [tau0, tau1.5] found in the best run. + +%Raghav Girgaonkar, April 2023 + +nSamples = length(inParams.dataX); + +fHandle = @(x) psofitfunc_tau_negative(x,inParams); + +params = inParams; + +nDim = 2; +outStruct = struct('bestLocation',[],... + 'bestFitness', [],... + 'totalFuncEvals',[],... + 'allBestFit',[],... + 'allBestLoc',[]); + +outResults = struct('allRunsOutput',struct('fitVal', [],... + 'qcCoefs',zeros(1,2),... + 'estTa',[],... + 'estSig',zeros(1,nSamples),... + 'totalFuncEvals',[],... + 'allBestFit',zeros(1,psoParams.maxSteps),... + 'allBestLoc',zeros(nDim,psoParams.maxSteps)),... + 'bestRun',[],... + 'bestFitness',[],... + 'bestSig', zeros(1,nSamples),... + 'bestQcCoefs',zeros(1,2),... + 'estAmp',[],... + 'estPhase',[],... + 'bestAmp',[],... + 'bestPhase',[],... + 'bestTime',[]); + +%Allocate storage for outputs: results from all runs are stored +for lpruns = 1:nRuns + outStruct(lpruns) = outStruct(1); + outResults.allRunsOutput(lpruns) = outResults.allRunsOutput(1); +end +%Independent runs of PSO in parallel. Change 'parfor' to 'for' if the +%parallel computing toolbox is not available. +% fprintf("Running PSO\n"); +parpool(nRuns); +parfor lpruns = 1:nRuns + %Reset random number generator for each worker + rng(lpruns); + outStruct(lpruns)=crcbpso(fHandle,nDim,psoParams,2); +end + +%Prepare output +fitVal = zeros(1,nRuns); +% sampling_interval = 1/sampling_freq; +%Time vectors for signal and total series +% timeVecSig = (0:(sampling_freq*T_sig - 1))/sampling_freq; +for lpruns = 1:nRuns + outResults.allRunsOutput(lpruns).allBestFit = outStruct(lpruns).allBestFit; + outResults.allRunsOutput(lpruns).allBestLoc = outStruct(lpruns).allBestLoc; + fitVal(lpruns) = outStruct(lpruns).bestFitness; + outResults.allRunsOutput(lpruns).fitVal = fitVal(lpruns); + [~,chirptimes,~] = fHandle(outStruct(lpruns).bestLocation); +% Q = qcCoefs +% index = ta_index + [estAmp, estTa, estPhase] = getparamestimates_negative(chirptimes, params); + outResults.allRunsOutput(lpruns).qcCoefs = chirptimes; + %Calculate time using sampling freq and ta_index +% estTa = ta_index/sampling_freq; + + outResults.allRunsOutput(lpruns).estTa = estTa; + tau0 = chirptimes(1); + tau1p5 = chirptimes(2); +% phaseq0 = gen2PNwaveform_tau_negative(params.fpos, estTa, 0, params.frange(1), params.frange(2), tau0,... +% tau1p5,params.datalen,0,1,params.N,params.avec, params.normfac); +% fftq0 = phaseq0; +% fftq1 = phaseq0.*params.phaseDiff; +% %Estimated Phase +% % yq0 = inParams.dataY*q0(:); +% % yq1 = inParams.dataY*q1(:); +% yq0 = innerprodpsd(fftq0, params.fftdataYbyPSD); +% yq1 = innerprodpsd(fftq1, params.fftdataYbyPSD); +% estPhase = atan2(yq1,yq0); + outResults.allRunsOutput(lpruns).estPhase = estPhase; + + %Estimated Amplitude +% estAmp = cos(estPhase)*yq0 + sin(estPhase)*yq1; + outResults.allRunsOutput(lpruns).estAmp = estAmp; + + %Estimated Signal +% estSigTemp = genqc(timeVecSig,1,qcCoefs,estPhase); + estSigphase = gen2PNwaveform_tau_negative(params.fpos, estTa, estPhase, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,estAmp,params.N,params.avec, params.normfac); + estSigfourier = (params.A).*estSigphase; + estSig = ifft(estSigfourier); +% estSigTemp_shifted = [zeros(1,floor(estTa*sampling_freq)-1), estSigTemp, zeros(1, nSamples - nSamplesSig - floor(estTa*sampling_freq)+1)]; +% estSig = estAmp*estSigTemp; + outResults.allRunsOutput(lpruns).estSig = estSig; + outResults.allRunsOutput(lpruns).totalFuncEvals = outStruct(lpruns).totalFuncEvals; +end +%Find the best run +[~,bestRun] = min(fitVal(:)); +outResults.bestRun = bestRun; +outResults.bestFitness = outResults.allRunsOutput(bestRun).fitVal; +outResults.bestSig = outResults.allRunsOutput(bestRun).estSig; +outResults.bestAmp = outResults.allRunsOutput(bestRun).estAmp; +outResults.bestPhase = outResults.allRunsOutput(bestRun).estPhase; +outResults.bestQcCoefs = outResults.allRunsOutput(bestRun).qcCoefs; +outResults.bestTime = outResults.allRunsOutput(bestRun).estTa; + diff --git a/Misc/NormalizationCheck/createPSD.m b/Misc/NormalizationCheck/createPSD.m new file mode 100644 index 0000000..60bd652 --- /dev/null +++ b/Misc/NormalizationCheck/createPSD.m @@ -0,0 +1,23 @@ +function [PSD, fvec]=createPSD(sampFreq, Tsig, welchPSD, freqs) +%% Script to create interpolated PSD from SHAPES and PWELCH Estimates + +%% Data Parameters +Fs = sampFreq; +T = Tsig; +N = Fs*T; + +kNyq = floor(N/2); +fvec = (0:(kNyq))*Fs/N; + +% nyqfreq = winlen*Fs; + +% freqs = (0:nyqfreq)*(Fs/(2*nyqfreq)); + +%% 1-D Interpolation +logwelchPSD = log10(welchPSD); +loginterPSD = interp1(freqs, logwelchPSD, fvec); + +% %% Antilog + +PSD = (10.^loginterPSD); +% PSD = (loginterPSD)./2; diff --git a/Misc/NormalizationCheck/files.json b/Misc/NormalizationCheck/files.json new file mode 100644 index 0000000..4b4908a --- /dev/null +++ b/Misc/NormalizationCheck/files.json @@ -0,0 +1,11 @@ +{ +"psdfile": "psdestimates.mat", +"noisefile": "50_noise_realizations.mat", +"datafile": "/Users/kqa493/Desktop/Accelerated-Network-Analysis/examples/L-L1_LOSC_CLN_4_V1-1187007040-2048_dtrndWhtnBndpss.hdf5", +"output_struct_location": "/scratch/09197/raghav/PSD-Estimation/outStruct_tau.mat", +"bestlocplot": "/scratch/09197/raghav/PSD-Estimation/bestloc_tau.png", +"bestfitplot": "/scratch/09197/raghav/PSD-Estimation/bestfitness_tau.png", +"psoresultplot": "/scratch/09197/raghav/PSD-Estimation/psoresults_tau.png", +"fitnessvaluesfile": "/scratch/09197/raghav/PSD-Estimation/fitness_values.txt", +"fitvalimgfile": "1-1.mat" +} diff --git a/Misc/NormalizationCheck/gen2PNtemplate_mass.m b/Misc/NormalizationCheck/gen2PNtemplate_mass.m new file mode 100644 index 0000000..ca6093b --- /dev/null +++ b/Misc/NormalizationCheck/gen2PNtemplate_mass.m @@ -0,0 +1,24 @@ +function signal = gen2PNtemplate_mass(fpos, ta, phase, fmin, fmax,m1,m2,datalen,initial_phase,snr,N, A, avec, PSDtotal) +%Creates and Returns normalized phase vector of waveform in the Fourier Domain +%General expression for a 2PN waveform in a fourier domain is +% W(f) = A(f)exp(-iPsi(f)), this function returns exp(-iPsi(f)) +% Waveform generation happens through chirp-time values tau0, tau1.5 + +fwavepos = waveform(fpos,ta,phase,fmin,fmax,m1,m2,datalen,initial_phase, avec); + +if mod(N,2) == 0 + fwaveneg = conj(fwavepos(end-1:-1:2)); +else + fwaveneg = conj(fwavepos(end:-1:2)); +end + + + +fwave = [fwavepos, fwaveneg]; +wf = ifft(A.*fwave); + +%Normalize to unit 1 +normfac = 1/sqrt(innerproduct(wf,wf,PSDtotal)); + +% Create final signal +signal = snr*normfac*wf; diff --git a/Misc/NormalizationCheck/gen2PNtemplate_tau.m b/Misc/NormalizationCheck/gen2PNtemplate_tau.m new file mode 100644 index 0000000..1cdaf96 --- /dev/null +++ b/Misc/NormalizationCheck/gen2PNtemplate_tau.m @@ -0,0 +1,27 @@ +function signal = gen2PNtemplate_tau(fpos, ta, phase, fmin, fmax,tau0, tau1p5,datalen,initial_phase,snr,N, A, avec, PSDtotal) +%Creates and Returns normalized phase vector of waveform in the Fourier Domain +%General expression for a 2PN waveform in a fourier domain is +% W(f) = A(f)exp(-iPsi(f)), this function returns exp(-iPsi(f)) +% Waveform generation happens through chirp-time values tau0, tau1.5 + +fwavepos = waveform_tau(fpos,ta,phase,fmin,fmax,tau0,tau1p5,datalen,initial_phase, avec); + +if mod(N,2) == 0 + fwaveneg = conj(fwavepos(end-1:-1:2)); +else + fwaveneg = conj(fwavepos(end:-1:2)); +end + + + +fwave = [fwavepos, fwaveneg]; +wf = ifft(A.*fwave); +fftwf = A.*fwave; + +%Normalize to unit 1 +% normfac = 1/sqrt(innerproduct_optmzd(fftwf,fftwf,PSDtotal)); +% normfac = 1/sqrt(innerproduct(wf,wf,PSDtotal)); +normfac = 1/sqrt(innerproduct_W(wf,wf,PSDtotal)); + +% Create final signal +signal = snr*normfac*wf; diff --git a/Misc/NormalizationCheck/gen2PNtemplate_tau_negative.m b/Misc/NormalizationCheck/gen2PNtemplate_tau_negative.m new file mode 100644 index 0000000..5086a55 --- /dev/null +++ b/Misc/NormalizationCheck/gen2PNtemplate_tau_negative.m @@ -0,0 +1,26 @@ +function signal = gen2PNtemplate_tau_negative(fpos, ta, phase, fmin, fmax,tau0, tau1p5,datalen,initial_phase,snr,N, A, avec, PSDtotal) +%Creates and Returns normalized phase vector of waveform in the Fourier Domain +%General expression for a 2PN waveform in a fourier domain is +% W(f) = A(f)exp(-iPsi(f)), this function returns exp(-iPsi(f)) +% Waveform generation happens through chirp-time values tau0, tau1.5 + +fwavepos = waveform_tau_negative(fpos,ta,phase,fmin,fmax,tau0,tau1p5,datalen,initial_phase, avec); + +if mod(N,2) == 0 + fwaveneg = conj(fwavepos(end-1:-1:2)); +else + fwaveneg = conj(fwavepos(end:-1:2)); +end + + + +fwave = [fwavepos, fwaveneg]; +wf = ifft(A.*fwave); +% fftwf = A.*fwave; + +%Normalize to unit 1 +% normfac = 1/sqrt(innerproduct_optmzd(fftwf,fftwf,PSDtotal)); +normfac = 1/sqrt(innerproduct_W(wf,wf,PSDtotal)); + +% Create final signal +signal = snr*normfac*wf; diff --git a/Misc/NormalizationCheck/gen2PNwaveform.m b/Misc/NormalizationCheck/gen2PNwaveform.m new file mode 100644 index 0000000..befa87f --- /dev/null +++ b/Misc/NormalizationCheck/gen2PNwaveform.m @@ -0,0 +1,20 @@ +function wave = gen2PNwaveform(fpos, ta, phase, fmin, fmax,m1, m2,datalen,initial_phase,snr,N,avec, normfac) +%Creates and Returns normalized phase vector of waveform in the Fourier Domain +%General expression for a 2PN waveform in a fourier domain is +% W(f) = A(f)exp(-iPsi(f)), this function returns exp(-iPsi(f)) +% Waveform generation happens through component mass values m1, m2 + + +fwavepos = waveform(fpos,ta,phase,fmin,fmax,m1,m2,datalen,initial_phase, avec); + +if mod(N,2) == 0 + fwaveneg = conj(fwavepos(end-1:-1:2)); +else + fwaveneg = conj(fwavepos(end:-1:2)); +end + +fwave = [fwavepos, fwaveneg]; + +wave = snr*normfac*fwave; + + diff --git a/Misc/NormalizationCheck/gen2PNwaveform_tau.m b/Misc/NormalizationCheck/gen2PNwaveform_tau.m new file mode 100644 index 0000000..8a31271 --- /dev/null +++ b/Misc/NormalizationCheck/gen2PNwaveform_tau.m @@ -0,0 +1,20 @@ +function wave = gen2PNwaveform_tau(fpos, ta, phase, fmin, fmax,tau0, tau1p5,datalen,initial_phase,snr,N, avec, normfac) +%Creates and Returns normalized phase vector of waveform in the Fourier Domain +%General expression for a 2PN waveform in a fourier domain is +% W(f) = A(f)exp(-iPsi(f)), this function returns exp(-iPsi(f)) +% Waveform generation happens through chirp-time values tau0, tau1.5 + +fwavepos = waveform_tau(fpos,ta,phase,fmin,fmax,tau0,tau1p5,datalen,initial_phase, avec); + +if mod(N,2) == 0 + fwaveneg = conj(fwavepos(end-1:-1:2)); +else + fwaveneg = conj(fwavepos(end:-1:2)); +end + +fwave = [fwavepos, fwaveneg]; + +wave = snr*normfac*fwave; + + + diff --git a/Misc/NormalizationCheck/gen2PNwaveform_tau_negative.m b/Misc/NormalizationCheck/gen2PNwaveform_tau_negative.m new file mode 100644 index 0000000..f3f54c4 --- /dev/null +++ b/Misc/NormalizationCheck/gen2PNwaveform_tau_negative.m @@ -0,0 +1,20 @@ +function wave = gen2PNwaveform_tau_negative(fpos, ta, phase, fmin, fmax,tau0, tau1p5,datalen,initial_phase,snr,N, avec, normfac) +%Creates and Returns normalized phase vector of waveform in the Fourier Domain +%General expression for a 2PN waveform in a fourier domain is +% W(f) = A(f)exp(-iPsi(f)), this function returns exp(-iPsi(f)) +% Waveform generation happens through chirp-time values tau0, tau1.5 + +fwavepos = waveform_tau_negative(fpos,ta,phase,fmin,fmax,tau0,tau1p5,datalen,initial_phase, avec); + +if mod(N,2) == 0 + fwaveneg = conj(fwavepos(end-1:-1:2)); +else + fwaveneg = conj(fwavepos(end:-1:2)); +end + +fwave = [fwavepos, fwaveneg]; + +wave = snr*normfac*fwave; + + + diff --git a/Misc/NormalizationCheck/gensignal.m b/Misc/NormalizationCheck/gensignal.m new file mode 100644 index 0000000..c5c0cda --- /dev/null +++ b/Misc/NormalizationCheck/gensignal.m @@ -0,0 +1,40 @@ +%% Script to generate a 2PN waveform normalized to a given SNR +function [signal]=gensignal() + %Injection parameters + masses = [4,300]; %Change to chirp-times as needed + initial_phase = 0; + phase = 0; %Coalescence phase of the signal + ta = 138; %TOA + datalen = 512; %Length of signal in seconds + frange = [30,700]; %Frequency range + snr = 15; % Injection SNR + sampFreq = 4096; %Sampling Frequency + N = datalen*sampFreq; %Total number of samples + fpos = (0:floor(datalen*sampFreq/2))*(1/datalen); %Positive Fourier frequency vector + + %PSD (Uncomment if using different PSD for normalization) + %Get two-sided PSD from iLIGO sensitivities + [~, PSD] = LIGOnoise(N,sampFreq,1,'sample'); + + %Make total PSD vector + negFStrt = 1-mod(N,2); + kNyq = floor(N/2)+1; + PSDtotal = [PSD, PSD((kNyq-negFStrt):-1:2)]; + + + %Preprocessing Vectors + [A,avec, phaseDiff] = preprocessing(frange(1), frange(2), fpos, datalen, datalen*sampFreq); + + + %Generate signal + %Comment below line and uncomment chirp time line if using chirptimes + signal = gen2PNtemplate_mass(fpos, ta, phase, frange(1), frange(2),masses(1), masses(2),datalen,initial_phase,snr,N, A, avec, PSDtotal); + + + %Uncomment if using chirptimes to create signal (Make sure to then comment + %above line) + % signal = gen2PNtemplate_tau(fpos, ta, phase, frange(1), frange(2),masses(1), masses(2),datalen,initial_phase,snr,N, A, avec, PSDtotal); + + %Uncomment if using negative chirp times to create waveform +% signal = gen2PNtemplate_tau_negative(fpos, ta, phase, frange(1), frange(2),masses(1), masses(2),datalen,initial_phase,snr,N, A, avec, PSDtotal); +end \ No newline at end of file diff --git a/Misc/NormalizationCheck/getparamestimates.m b/Misc/NormalizationCheck/getparamestimates.m new file mode 100644 index 0000000..5e70856 --- /dev/null +++ b/Misc/NormalizationCheck/getparamestimates.m @@ -0,0 +1,31 @@ +function [estAmp, estTa, estPhase] = getparamestimates(x, params) +%UNTITLED Summary of this function goes here +% Detailed explanation goes here + +%Generate normalized quadrature templates +tau0 = x(1); +tau1p5 = x(2); +fftq0 = gen2PNtemplate_tau(params.fpos, 0, 0, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.A,params.avec, params.PSDtotal); +fftq1 = gen2PNtemplate_tau(params.fpos, 0, pi/2, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.A,params.avec, params.PSDtotal); + +%Compute fitness value after maximizing by matched filtering +% mf1 = mlftr(params.dataY, q0, params.PSDtotal); +% mf2 = mlftr(params.dataY, q1, params.PSDtotal); +% mf1 = mlftr2(params.dataY, q0, params.TFtotal); +% mf2 = mlftr2(params.dataY, q1, params.TFtotal); +mf1 = matchedfiltering(params.fftdataYbyPSD, fftq0); +mf2 = matchedfiltering(params.fftdataYbyPSD, fftq1); +[max_val, max_arg] = max(mf1(1:end - params.T_sig*params.Fs).^2 + mf2(1:end - params.T_sig*params.Fs).^2); +%Estimated SNR +estAmp = sqrt(max_val); +%Estimated TOA: +estTa = max_arg/params.Fs; + +yq0 = mf1(max_arg); +yq1 = mf2(max_arg); + +estPhase = atan2(yq1,yq0); +end + diff --git a/Misc/NormalizationCheck/getparamestimates_mass.m b/Misc/NormalizationCheck/getparamestimates_mass.m new file mode 100644 index 0000000..c327310 --- /dev/null +++ b/Misc/NormalizationCheck/getparamestimates_mass.m @@ -0,0 +1,26 @@ +function [estAmp, estTa, estPhase] = getparamestimates_mass(masses, params) + + m1 = masses(1); + m2 = masses(2); + phaseq0 = gen2PNwaveform(params.fpos, 0, 0, params.frange(1), params.frange(2), m1,... + m2,params.datalen,0,1,params.N,params.avec, params.normfac); + + fftq0 = phaseq0; + + fftq1 = phaseq0.*params.phaseDiff; + + %Compute fitness value after maximizing by matched filtering + mf1 = matchedfiltering(params.fftdataYbyPSD, fftq0); + mf2 = matchedfiltering(params.fftdataYbyPSD, fftq1); + [max_val, max_arg] = max(mf1(1:end - params.T_sig*params.Fs).^2 + mf2(1:end - params.T_sig*params.Fs).^2); + %Estimated SNR + estAmp = sqrt(max_val); + %Estimated TOA: + estTa = max_arg/params.Fs; + + yq0 = mf1(max_arg); + yq1 = mf2(max_arg); + + estPhase = atan2(yq1,yq0); +end + diff --git a/Misc/NormalizationCheck/getparamestimates_negative.m b/Misc/NormalizationCheck/getparamestimates_negative.m new file mode 100644 index 0000000..0bfb5f9 --- /dev/null +++ b/Misc/NormalizationCheck/getparamestimates_negative.m @@ -0,0 +1,28 @@ +function [estAmp, estTa, estPhase] = getparamestimates_negative(chirptimes, params) +%UNTITLED Summary of this function goes here +% Detailed explanation goes here + +tau0 = chirptimes(1); +tau1p5 = chirptimes(2); +phaseq0 = gen2PNwaveform_tau_negative(params.fpos, 0, 0, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.avec, params.normfac); + +fftq0 = phaseq0; + +fftq1 = phaseq0.*params.phaseDiff; + +%Compute fitness value after maximizing by matched filtering +mf1 = matchedfiltering(params.fftdataYbyPSD, fftq0); +mf2 = matchedfiltering(params.fftdataYbyPSD, fftq1); +[max_val, max_arg] = max(mf1(1:end - params.T_sig*params.Fs).^2 + mf2(1:end - params.T_sig*params.Fs).^2); +%Estimated SNR +estAmp = sqrt(max_val); +%Estimated TOA: +estTa = max_arg/params.Fs; + +yq0 = mf1(max_arg); +yq1 = mf2(max_arg); + +estPhase = atan2(yq1,yq0); +end + diff --git a/Misc/NormalizationCheck/iLIGOSensitivity.txt b/Misc/NormalizationCheck/iLIGOSensitivity.txt new file mode 100755 index 0000000..3fba188 --- /dev/null +++ b/Misc/NormalizationCheck/iLIGOSensitivity.txt @@ -0,0 +1 @@ +1 1.57837977 1.1 0.415636725 1.21 0.10945014 1.33 0.029126516 1.46 0.007893535 1.61 0.002007471 1.77 0.000532826 1.95 0.000137318 2.14 3.74E-05 2.36 9.49E-06 2.59 2.58E-06 2.85 6.77E-07 3.14 1.74E-07 3.45 4.66E-08 3.8 1.21E-08 4.18 3.17E-09 4.59 8.57E-10 5.05 2.25E-10 5.56 5.85E-11 6.12 1.53E-11 6.73 4.04E-12 7.4 1.07E-12 8.14 2.81E-13 8.95 7.46E-14 9.85 1.95E-14 10.83 5.17E-15 11.92 1.35E-15 13.11 3.56E-16 14.42 9.39E-17 15.86 2.48E-17 17.45 6.50E-18 19.19 1.72E-18 21.11 4.52E-19 23.23 1.18E-19 25.55 3.12E-20 28.1 8.26E-21 30.91 2.21E-21 34 6.60E-22 37.4 3.11E-22 41.14 2.28E-22 45.26 1.86E-22 49.79 1.54E-22 54.76 1.27E-22 60.24 1.06E-22 66.26 8.77E-23 72.89 7.30E-23 80.18 6.10E-23 88.2 5.13E-23 97.02 4.37E-23 106.72 3.77E-23 117.39 3.32E-23 129.13 3.02E-23 142.04 2.83E-23 156.25 2.75E-23 171.87 2.75E-23 189.06 2.84E-23 207.97 2.98E-23 228.76 3.17E-23 251.64 3.41E-23 276.8 3.69E-23 304.48 4.01E-23 334.93 4.37E-23 368.42 4.77E-23 405.27 5.22E-23 445.79 5.71E-23 490.37 6.26E-23 539.41 6.87E-23 593.35 7.54E-23 652.68 8.27E-23 717.95 9.09E-23 789.75 9.98E-23 868.72 1.10E-22 955.59 1.21E-22 1051.15 1.32E-22 1156.27 1.46E-22 1271.9 1.60E-22 1399.08 1.76E-22 1538.99 1.94E-22 1692.89 2.13E-22 1862.18 2.34E-22 2048.4 2.57E-22 2253.24 2.83E-22 2478.56 3.11E-22 2726.42 3.43E-22 2999.06 3.77E-22 3298.97 4.14E-22 3628.87 4.56E-22 3991.75 5.01E-22 4390.93 5.51E-22 4830.02 6.07E-22 5313.02 6.67E-22 5844.32 7.34E-22 6428.76 8.07E-22 7071.63 8.88E-22 7778.8 9.77E-22 8556.68 1.07E-21 9412.34 1.18E-21 \ No newline at end of file diff --git a/Misc/NormalizationCheck/injsnr_histogram.png b/Misc/NormalizationCheck/injsnr_histogram.png new file mode 100644 index 0000000..603cf3b Binary files /dev/null and b/Misc/NormalizationCheck/injsnr_histogram.png differ diff --git a/Misc/NormalizationCheck/innerprodpsd.m b/Misc/NormalizationCheck/innerprodpsd.m new file mode 100755 index 0000000..a499f0a --- /dev/null +++ b/Misc/NormalizationCheck/innerprodpsd.m @@ -0,0 +1,15 @@ +function innProd = innerprodpsd(fftY,fftXbyPSD) +%Calculates the inner product of vectors X and Y for the case of Gaussian +% stationary noise having a specified power spectral density (PSD) Sn(f). +% Input: fftY = Normlized FFT vector of the phase term of quadrature template +% fftXbyPSD = FFT vector of data divided by PSD Sn(f) and +% multiplied by the frequency amplitude term A. +% Output: innProd: inner product value +%The scalar normalization factor (1/N) is derived from Parseval's theorem. + +%Raghav Girgaonkar, Apr 2023 + +nSamples = length(fftXbyPSD); +dataLen = nSamples; +innProd = (1/dataLen)*(fftXbyPSD)*fftY'; +innProd = real(innProd); diff --git a/Misc/NormalizationCheck/innerproduct.m b/Misc/NormalizationCheck/innerproduct.m new file mode 100644 index 0000000..289c1f0 --- /dev/null +++ b/Misc/NormalizationCheck/innerproduct.m @@ -0,0 +1,10 @@ +function [innerProd] = innerproduct(X,Y,PSD, sampFreq) +%INNERPRODUCT +% N = length(X); +N = sampFreq*length(X); +% Fs = 4096; +% N = length(X)/Fs; +innerProd = (1/N)*sum((fft(X)./PSD).*conj(fft(Y))); +innerProd = real(innerProd); +end + diff --git a/Misc/NormalizationCheck/innerproduct_W.m b/Misc/NormalizationCheck/innerproduct_W.m new file mode 100644 index 0000000..1178c88 --- /dev/null +++ b/Misc/NormalizationCheck/innerproduct_W.m @@ -0,0 +1,10 @@ +function [innerProd] = innerproduct_W(X,Y,PSD) +%INNERPRODUCT + N = length(X); +% N = 4096*length(X); +% Fs = 4096; +% N = length(X)/Fs; +innerProd = (1/N)*sum((fft(X)./PSD).*conj(fft(Y))); +innerProd = real(innerProd); +end + diff --git a/Misc/NormalizationCheck/innerproduct_optmzd.m b/Misc/NormalizationCheck/innerproduct_optmzd.m new file mode 100644 index 0000000..c4f8a1f --- /dev/null +++ b/Misc/NormalizationCheck/innerproduct_optmzd.m @@ -0,0 +1,7 @@ +function [innerProd] = innerproduct_optmzd(X,Y,PSD) +%INNERPRODUCT +N = length(X); +innerProd = (1/N)*sum((X./PSD).*conj(Y)); + +end + diff --git a/Misc/NormalizationCheck/matchedfiltering.m b/Misc/NormalizationCheck/matchedfiltering.m new file mode 100644 index 0000000..7e284d0 --- /dev/null +++ b/Misc/NormalizationCheck/matchedfiltering.m @@ -0,0 +1,12 @@ +function timesVec = matchedfiltering(fftXbyPSD, fftY) +%Function to perform FFT correlation based matched filtering +% Input: fftXbyPSD = FFT vector of data which has been multiplied by +% frequency magnitude vector A and divided by total PSD +% fftY = Normalized FFT phase vector of quadrature template +% Output: timesVec = MatchedFiltering timeseries vector + +%Raghav Girgaonkar, April 2023 + +innProd = fftXbyPSD.*conj(fftY); +timesVec = ifft(innProd); + diff --git a/Misc/NormalizationCheck/mfgw_mass.m b/Misc/NormalizationCheck/mfgw_mass.m new file mode 100644 index 0000000..440def2 --- /dev/null +++ b/Misc/NormalizationCheck/mfgw_mass.m @@ -0,0 +1,27 @@ +function [mfVal, max_arg] = mfgw_mass(x,params) +%MatchedFiltering for Mass space PSO +%Generates a combined matched-filtering timeseries from both quadrature templates and +%returns index and value of the maximum of this series, +%Input: x = [m1, m2], vector containing mass parameters for +% creating quadrature templates +% params: Struct containing signal parameters +%Output: mfVal: Maximum value of total matchedfiltering timeseries +% max_arg: Index of maximum value + + +%Raghav Girgaonkar, April 2023 +m1 = x(1); +m2 = x(2); +phaseq0 = gen2PNwaveform(params.fpos, 0, 0, params.frange(1), params.frange(2), m1,... + m2,params.datalen,0,1,params.N,params.avec, params.normfac); + +fftq0 = phaseq0; + +fftq1 = phaseq0.*params.phaseDiff; + + +%Compute fitness value after maximizing by matched filtering +mf1 = matchedfiltering(params.fftdataYbyPSD, fftq0); +mf2 = matchedfiltering(params.fftdataYbyPSD, fftq1); +[max_val, max_arg] = max(mf1(1:end - params.T_sig*params.Fs).^2 + mf2(1:end - params.T_sig*params.Fs).^2); +mfVal = -1*max_val; diff --git a/Misc/NormalizationCheck/mfgw_tau.m b/Misc/NormalizationCheck/mfgw_tau.m new file mode 100644 index 0000000..3ddc98c --- /dev/null +++ b/Misc/NormalizationCheck/mfgw_tau.m @@ -0,0 +1,27 @@ +function [mfVal, max_arg] = mfgw_tau(x,params) +%MatchedFiltering for Chirp time space PSO +%Generates a combined matched-filtering timeseries from both quadrature templates and +%returns index and value of the maximum of this series, +%Input: x = [tau0, tau1.5], vector containing chirp-time parameters for +% creating quadrature templates +% params: Struct containing signal parameters +%Output: mfVal: Maximum value of total matchedfiltering timeseries +% max_arg: Index of maximum value + +%Raghav Girgaonkar, April 2023 + +%Generate normalized quadrature templates +tau0 = x(1); +tau1p5 = x(2); +phaseq0 = gen2PNwaveform_tau(params.fpos, 0, 0, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.avec, params.normfac); + +fftq0 = phaseq0; + +fftq1 = phaseq0.*params.phaseDiff; + +%Compute fitness value after maximizing by matched filtering +mf1 = matchedfiltering(params.fftdataYbyPSD, fftq0); +mf2 = matchedfiltering(params.fftdataYbyPSD, fftq1); +[max_val, max_arg] = max(mf1(1:end - params.T_sig*params.Fs).^2 + mf2(1:end - params.T_sig*params.Fs).^2); +mfVal = -1*max_val; diff --git a/Misc/NormalizationCheck/mfgw_tau_negative.m b/Misc/NormalizationCheck/mfgw_tau_negative.m new file mode 100644 index 0000000..ffeaf06 --- /dev/null +++ b/Misc/NormalizationCheck/mfgw_tau_negative.m @@ -0,0 +1,27 @@ +function [mfVal, max_arg] = mfgw_tau_negative(x,params) +%MatchedFiltering for Chirp time space PSO +%Generates a combined matched-filtering timeseries from both quadrature templates and +%returns index and value of the maximum of this series, +%Input: x = [tau0, tau1.5], vector containing chirp-time parameters for +% creating quadrature templates +% params: Struct containing signal parameters +%Output: mfVal: Maximum value of total matchedfiltering timeseries +% max_arg: Index of maximum value + +%Raghav Girgaonkar, April 2023 + +%Generate normalized quadrature templates +tau0 = x(1); +tau1p5 = x(2); +phaseq0 = gen2PNwaveform_tau_negative(params.fpos, 0, 0, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.avec, params.normfac); + +fftq0 = phaseq0; + +fftq1 = phaseq0.*params.phaseDiff; + +%Compute fitness value after maximizing by matched filtering +mf1 = matchedfiltering(params.fftdataYbyPSD, fftq0); +mf2 = matchedfiltering(params.fftdataYbyPSD, fftq1); +[max_val, max_arg] = max(mf1(1:end - params.T_sig*params.Fs).^2 + mf2(1:end - params.T_sig*params.Fs).^2); +mfVal = -1*max_val; diff --git a/Misc/NormalizationCheck/mfgw_tau_new.m b/Misc/NormalizationCheck/mfgw_tau_new.m new file mode 100644 index 0000000..6bdbddb --- /dev/null +++ b/Misc/NormalizationCheck/mfgw_tau_new.m @@ -0,0 +1,28 @@ +function [mf1, mf2, mfVal, max_arg] = mfgw_tau_new(x,params) +%MatchedFiltering for Chirp time space PSO +%Generates a combined matched-filtering timeseries from both quadrature templates and +%returns index and value of the maximum of this series, +%Input: x = [tau0, tau1.5], vector containing chirp-time parameters for +% creating quadrature templates +% params: Struct containing signal parameters +%Output: mfVal: Maximum value of total matchedfiltering timeseries +% max_arg: Index of maximum value + +%Raghav Girgaonkar, April 2023 + +%Generate normalized quadrature templates +tau0 = x(1); +tau1p5 = x(2); +q0 = gen2PNtemplate_tau(params.fpos, 0, 0, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.A,params.avec, params.PSDtotal); +q1 = gen2PNtemplate_tau(params.fpos, 0, pi/2, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.A,params.avec, params.PSDtotal); + +%Compute fitness value after maximizing by matched filtering +mf1 = mlftr(params.dataY, q0, params.PSDtotal); +mf2 = mlftr(params.dataY, q1, params.PSDtotal); +% mf1 = mlftr2(params.dataY, q0, params.TFtotal); +% mf2 = mlftr2(params.dataY, q1, params.TFtotal); +% [max_val, max_arg] = max(mf1(1:end - params.T_sig*params.Fs).^2 + mf2(1:end - params.T_sig*params.Fs).^2); +[max_val, max_arg] = max(mf1(1:end - params.T_sig*params.Fs).^2 + mf2(1:end - params.T_sig*params.Fs).^2); +mfVal = -1*max_val; diff --git a/Misc/NormalizationCheck/mfgw_tau_whtnd.m b/Misc/NormalizationCheck/mfgw_tau_whtnd.m new file mode 100644 index 0000000..53c089f --- /dev/null +++ b/Misc/NormalizationCheck/mfgw_tau_whtnd.m @@ -0,0 +1,27 @@ +function [mfVal, max_arg] = mfgw_tau_whtnd(x,params) +%MatchedFiltering for Chirp time space PSO +%Generates a combined matched-filtering timeseries from both quadrature templates and +%returns index and value of the maximum of this series, +%Input: x = [tau0, tau1.5], vector containing chirp-time parameters for +% creating quadrature templates +% params: Struct containing signal parameters +%Output: mfVal: Maximum value of total matchedfiltering timeseries +% max_arg: Index of maximum value + +%Raghav Girgaonkar, April 2023 + +%Generate normalized quadrature templates +tau0 = x(1); +tau1p5 = x(2); +fftq0 = gen2PNtemplate_tau(params.fpos, 0, 0, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.A,params.avec, params.PSDtotal); +fftq1 = gen2PNtemplate_tau(params.fpos, 0, pi/2, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.A,params.avec, params.PSDtotal); + +%Compute fitness value after maximizing by matched filtering +% mf1 = mlftr(params.dataY, q0, params.PSDtotal); +% mf2 = mlftr(params.dataY, q1, params.PSDtotal); +mf1 = matchedfiltering(params.fftdataYbyPSD, fftq0); +mf2 = matchedfiltering(params.fftdataYbyPSD, fftq1); +[max_val, max_arg] = max(mf1(1:end - params.T_sig*params.Fs).^2 + mf2(1:end - params.T_sig*params.Fs).^2); +mfVal = -1*max_val; diff --git a/Misc/NormalizationCheck/mfgw_tau_whtnd_og.m b/Misc/NormalizationCheck/mfgw_tau_whtnd_og.m new file mode 100644 index 0000000..b4ca9fa --- /dev/null +++ b/Misc/NormalizationCheck/mfgw_tau_whtnd_og.m @@ -0,0 +1,27 @@ +function [mf1,mf2,mfVal, max_arg] = mfgw_tau_whtnd_og(x,params) +%MatchedFiltering for Chirp time space PSO +%Generates a combined matched-filtering timeseries from both quadrature templates and +%returns index and value of the maximum of this series, +%Input: x = [tau0, tau1.5], vector containing chirp-time parameters for +% creating quadrature templates +% params: Struct containing signal parameters +%Output: mfVal: Maximum value of total matchedfiltering timeseries +% max_arg: Index of maximum value + +%Raghav Girgaonkar, April 2023 + +%Generate normalized quadrature templates +tau0 = x(1); +tau1p5 = x(2); +q0 = gen2PNtemplate_tau(params.fpos, 0, 0, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.A,params.avec, params.PSDtotal); +q1 = gen2PNtemplate_tau(params.fpos, 0, pi/2, params.frange(1), params.frange(2), tau0,... + tau1p5,params.datalen,0,1,params.N,params.A,params.avec, params.PSDtotal); + +%Compute fitness value after maximizing by matched filtering +mf1 = mlftr2(params.dataY, q0, params.TFtotal); +mf2 = mlftr2(params.dataY, q1, params.TFtotal); +% mf1 = matchedfiltering(params.fftdataYbyPSD, fftq0); +% mf2 = matchedfiltering(params.fftdataYbyPSD, fftq1); +[max_val, max_arg] = max(mf1(1:end - params.T_sig*params.Fs).^2 + mf2(1:end - params.T_sig*params.Fs).^2); +mfVal = -1*max_val; diff --git a/Misc/NormalizationCheck/mlftr.m b/Misc/NormalizationCheck/mlftr.m new file mode 100644 index 0000000..307b64a --- /dev/null +++ b/Misc/NormalizationCheck/mlftr.m @@ -0,0 +1,5 @@ +function [mftimeseries] = mlftr(X,Y,PSD) +%MLFTR +mftimeseries = (1/4096)*ifft((fft(X)./PSD).*conj(fft(Y))); +end + diff --git a/Misc/NormalizationCheck/mlftr2.m b/Misc/NormalizationCheck/mlftr2.m new file mode 100644 index 0000000..c17930b --- /dev/null +++ b/Misc/NormalizationCheck/mlftr2.m @@ -0,0 +1,5 @@ +function [mftimeseries] = mlftr2(X,Y,TFtotal) +%MLFTR on whitened X +mftimeseries = ifft((fft(X).*TFtotal).*conj(fft(Y))); +end + diff --git a/Misc/NormalizationCheck/optimizedMatchedFiltering.m b/Misc/NormalizationCheck/optimizedMatchedFiltering.m new file mode 100644 index 0000000..1092f30 --- /dev/null +++ b/Misc/NormalizationCheck/optimizedMatchedFiltering.m @@ -0,0 +1,174 @@ +%% Script to demonstrate code optimization + +%% Define Data Parameters +sampFreq = 4096; +datalen = 512; +N = datalen*sampFreq; +frange = [30,700]; +fpos = (0:floor(N/2))*(1/datalen); +timeVec = (0:(datalen*sampFreq - 1))/sampFreq; + +%% Create Colored Gaussian Noise +[noise, PSD] = LIGOnoise(N, sampFreq, 1, 'sample'); +% S = load('realization_512s_inj138_fs4096_n1.mat'); +% noise = S.dataY; +% PSD = S.dsstPSD; + +%Create Total PSD and Transfer Function +negFStrt = 1-mod(N,2); +kNyq = floor(N/2)+1; + +PSDtotal = [PSD, PSD((kNyq-negFStrt):-1:2)]; +TFtotal = 1./sqrt(PSDtotal); + + +%Pre-calculate frequency magnitude and phase difference terms + [A,avec, phaseDiff] = preprocessing(frange(1),frange(2),fpos, datalen, N); + +%Create general normalization factor N_fW +AbysqrtPSD = A.*TFtotal; +innProd = (1/N)*(AbysqrtPSD)*AbysqrtPSD'; +genNormfacSqr = real(innProd); +N_fW = 1/sqrt(genNormfacSqr); + +%% Generate Signal to Inject + chirptimes = [29.6373, 1.1045]; + ta = 138; + initial_phase = 0; + phase = 0; + snr = 30; + + + + % Generate unnormalized signal vector of arbitrary SNR +fwavepos = waveform_tau(fpos,ta,phase,frange(1),frange(2),chirptimes(1),chirptimes(2),datalen,initial_phase, avec); + +if mod(N,2) == 0 + fwaveneg = conj(fwavepos(end-1:-1:2)); +else + fwaveneg = conj(fwavepos(end:-1:2)); +end + +fwave = [fwavepos, fwaveneg]; +signal = ifft(A.*fwave); + + + +%Normalize signal to preset SNR using expression for strain innerproduct +normFac = real(snr*(1/sqrt((1/(N*sampFreq))*(sum((fft(signal).*TFtotal).*conj(fft(signal).*TFtotal)))))); +innProd = (1/N)*(fft(signal).*TFtotal)*(fft(signal).*TFtotal)'; +normFacW = snr/sqrt(innProd); +% normFacW = real(snr*(1/sqrt((1/(N))*(sum((fft(signal).*TFtotal).*conj(fft(signal).*TFtotal)))))); + + +signal = signal*normFac; +% signal = signal*normFacW; +% signal = signal*N_fW; + +% whtndsignal = ifft(fft(signal).*TFtotal); + +% strain = noise; + +%Inject signal directly into strain + strain = noise + signal; +% strain = filtdata + signal; + + +%Highpassing (Optional) +% rolloff = 0.125; +% fmin = frange(1); +% segwin = strain.*tukeywin(length(strain),rolloff*sampFreq/length(strain))'; +% seghpass = highpass(segwin, fmin, sampFreq, ImpulseResponse="iir",Steepness=0.95); +% strain = seghpass; + +%Tukey Window data +rolloff = 0.5; %Roll-off in seconds +winstrain = strain.*tukeywin(length(strain), rolloff*sampFreq/N)'; + +% mf = (1/sampFreq)*ifft((fft(winstrain).*TFtotal).*conj(fft(qsignal).*TFtotal)); + + + +%% Whiten and Normalize Strain data + +whtndstrain = ifft((1/sqrt(sampFreq))*fft(winstrain).*TFtotal); + +% whtndstrain = whtndstrain + whtndsignal; + +%% Perform Optimized Matched Filtering + +%Create fftdatabyPSD vector +fftdatabyPSD = (A.*fft(whtndstrain).*TFtotal); + + +% N_fW = 1/sqrt((1/N)*sum((A.*TFtotal).*conj(A.*TFtotal))); + +%Generate phase Fourier terms of quadrature templates q0,q1 using the same +%parameters used for signal injection +fwavepos = waveform_tau(fpos,0,phase,frange(1),frange(2),chirptimes(1),chirptimes(2),datalen,initial_phase, avec); + +if mod(N,2) == 0 + fwaveneg = conj(fwavepos(end-1:-1:2)); +else + fwaveneg = conj(fwavepos(end:-1:2)); +end + +phaseq0 = [fwavepos, fwaveneg]; + +phaseq1 = phaseq0.*phaseDiff; %Phase difference of pi/2 + +%Calculate M1 and M2 +m1 = ifft(fftdatabyPSD.*conj(N_fW*phaseq0)); +m2 = ifft(fftdatabyPSD.*conj(N_fW*phaseq1)); + +%Generate final timeseries +m = sqrt(m1.^2 + m2.^2); + +%Calculate estimate Time of Arrival and SNR +[SNR,tindex] = max(m); +TOA = tindex/sampFreq; +% SNRs(i) = SNR; +% TOAs(i) = TOA; + + +% figure; +% histogram(SNRs,'NumBins',50); +% +% +% figure; +% histogram(TOAs,'NumBins',50); + + +% %Display parameters +dispstr = ['The estimated TOA = ',num2str(TOA),' and the estimated SNR = ',num2str(SNR)]; +disp(dispstr); + +%% Plots +%Plot m1 and m2 individually +figure; +plot(timeVec, m1); +xlabel('Time (s)'); +title('M_1 Matched Filter Timeseries'); + +figure; +plot(timeVec, m2); +xlabel('Time (s)'); +title('M_2 Matched Filter Timeseries'); + +figure; +plot(timeVec, m); +xlabel('Time (s)'); +title('$\sqrt{M_1^2 + M_2^2}$ Matched Filter Timeseries','Interpreter','latex'); + + + + + + + + + + + + + diff --git a/Misc/NormalizationCheck/preprocessing.m b/Misc/NormalizationCheck/preprocessing.m new file mode 100644 index 0000000..244b63a --- /dev/null +++ b/Misc/NormalizationCheck/preprocessing.m @@ -0,0 +1,46 @@ +function [A,avec, phaseDiff] = preprocessing(fmin,fmax,fpos, datalen, N) +%PREPROCESSING Function to create vectors used for waveform generation +% Input: fmin: Low frequency cutoff for waveform +% fmax: High frequency cutoff for waveform +% fpos: Positive frequency vector +% datalen: total length of data in seconds +% N: Total number of samples in data segment +% Output: A: Frequency magnitude vector +% avec: matrix containing alpha terms used for waveform +% generation +% phaseDiff: phase difference vector for quadrature templates + +%Raghav Girgaonkar, May 2023 + + +%Make Frequency Magnitude and Phase difference vectors +Apos = zeros(size(fpos)); +phaseDiffpos = -1j*ones(size(fpos)); + +Apos(2:end) = fpos(2:end).^(-7/6); +%Modify Apos +min_index = floor(datalen*fmin) + 1; +max_index = floor(datalen*fmax) + 1; +Apos(1:min_index-1) = 0; +Apos(max_index + 1: end) = 0; +%Make Aneg and Phasediffneg +if mod(N,2) == 0 + Aneg = conj(Apos(end-1:-1:2)); + phaseDiffneg = conj(phaseDiffpos(end-1:-1:2)); +else + Aneg = conj(Apos(end:-1:2)); + phaseDiffneg = conj(phaseDiffpos(end:-1:2)); +end +%Make full A +A = [Apos, Aneg]; +phaseDiff = [phaseDiffpos, phaseDiffneg]; + +a0fvec = ((fpos(2:end)./fmin).^(-5/3)); +a1fvec = ((fpos(2:end)./fmin).^(-4/3)); +a2fvec = ((fpos(2:end)./fmin).^(-3/3)); +a3fvec = ((fpos(2:end)./fmin).^(-2/3)); +a4fvec = ((fpos(2:end)./fmin).^(-1/3)); + +avec = [a0fvec ; a1fvec; a2fvec; a3fvec; a4fvec]; +end + diff --git a/Misc/NormalizationCheck/pso.json b/Misc/NormalizationCheck/pso.json new file mode 100644 index 0000000..6b0a5f1 --- /dev/null +++ b/Misc/NormalizationCheck/pso.json @@ -0,0 +1,5 @@ +{ +"type": "tau", +"maxSteps": 1, +"nruns": 1 +} diff --git a/Misc/NormalizationCheck/psofitfunc.m b/Misc/NormalizationCheck/psofitfunc.m new file mode 100644 index 0000000..dc6cd5d --- /dev/null +++ b/Misc/NormalizationCheck/psofitfunc.m @@ -0,0 +1,39 @@ +function [fitVal,varargout] = psofitfunc(xVec,params) +%Fitness function for Mass Space PSO +%Input: xVec: normalized location vector of mass parameters +% params: Struct containing signal and data parameters +%Output: fitval: Fitness value at location specified by xVec +% varargout: Additional output arguments sucha as the index of max +% value in matchedfiltering timeseries + +%Raghav Girgaonkar, April 2023 + +%rows: points +%columns: coordinates of a point +[nVecs,~]=size(xVec); + +%storage for fitness values +fitVal = zeros(nVecs,1); + +%Check for out of bound coordinates and flag them +validPts = crcbchkstdsrchrng(xVec); +%Set fitness for invalid points to infty +fitVal(~validPts)=inf; +xVec(validPts,:) = s2rv(xVec(validPts,:),params); + +for lpc = 1:nVecs + if validPts(lpc) + % Only the body of this block should be replaced for different fitness + % functions + x = xVec(lpc,:); + [fitVal(lpc), max_index] = mfgw_mass(x, params); + end +end + +%Return max_index if requested +if nargout > 1 + varargout{1} = xVec; + varargout{2}=max_index; +end + + diff --git a/Misc/NormalizationCheck/psofitfunc_tau.m b/Misc/NormalizationCheck/psofitfunc_tau.m new file mode 100644 index 0000000..75a7e51 --- /dev/null +++ b/Misc/NormalizationCheck/psofitfunc_tau.m @@ -0,0 +1,38 @@ +function [fitVal,varargout] = psofitfunc_tau(xVec,params) +%Fitness function for Chirp-time Space PSO +%Input: xVec: normalized location vector of mass parameters +% params: Struct containing signal and data parameters +%Output: fitval: Fitness value at location specified by xVec +% varargout: Additional output arguments sucha as the index of max +% value in matchedfiltering timeseries + +%Raghav Girgaonkar, April 2023 + +%rows: points +%columns: coordinates of a point +[nVecs,~]=size(xVec); + +%storage for fitness values +fitVal = zeros(nVecs,1); + +%Check for out of bound coordinates and flag them +validPts = crcbchkstdsrchrng(xVec); +%Set fitness for invalid points to infty +fitVal(~validPts)=inf; +xVec(validPts,:) = s2rv(xVec(validPts,:),params); + +for lpc = 1:nVecs + if validPts(lpc) + % Only the body of this block should be replaced for different fitness + % functions + x = xVec(lpc,:); + [fitVal(lpc), max_index] = mfgw_tau_whtnd(x, params); + end +end + +%Return max_index if requested +if nargout > 1 + varargout{1} = xVec; + varargout{2}=max_index; +end + diff --git a/Misc/NormalizationCheck/psofitfunc_tau_negative.m b/Misc/NormalizationCheck/psofitfunc_tau_negative.m new file mode 100644 index 0000000..3074f6e --- /dev/null +++ b/Misc/NormalizationCheck/psofitfunc_tau_negative.m @@ -0,0 +1,38 @@ +function [fitVal,varargout] = psofitfunc_tau_negative(xVec,params) +%Fitness function for Chirp-time Space PSO +%Input: xVec: normalized location vector of mass parameters +% params: Struct containing signal and data parameters +%Output: fitval: Fitness value at location specified by xVec +% varargout: Additional output arguments sucha as the index of max +% value in matchedfiltering timeseries + +%Raghav Girgaonkar, April 2023 + +%rows: points +%columns: coordinates of a point +[nVecs,~]=size(xVec); + +%storage for fitness values +fitVal = zeros(nVecs,1); + +%Check for out of bound coordinates and flag them +validPts = crcbchkstdsrchrng(xVec); +%Set fitness for invalid points to infty +fitVal(~validPts)=inf; +xVec(validPts,:) = s2rv(xVec(validPts,:),params); + +for lpc = 1:nVecs + if validPts(lpc) + % Only the body of this block should be replaced for different fitness + % functions + x = xVec(lpc,:); + [fitVal(lpc), max_index] = mfgw_tau_negative(x, params); + end +end + +%Return max_index if requested +if nargout > 1 + varargout{1} = xVec; + varargout{2}=max_index; +end + diff --git a/Misc/NormalizationCheck/rungwpso.m b/Misc/NormalizationCheck/rungwpso.m new file mode 100644 index 0000000..b96497b --- /dev/null +++ b/Misc/NormalizationCheck/rungwpso.m @@ -0,0 +1,326 @@ +function []=rungwpso(filename) +%Script to run Chirp-time or Mass Space PSO on a user specified datafile and PSD. +%Refer to User Manual for usage + +%Raghav Girgaonkar, Apr 2023 + +%% Read JSON Files +fname = filename; +str = fileread(fname); +filenames = jsondecode(str); +fname = filenames.signalparamfile; +str = fileread(fname); +params = jsondecode(str); +fname = filenames.psoparamfile; +str = fileread(fname); +pso = jsondecode(str); +fname = filenames.filenames; +str = fileread(fname); +files = jsondecode(str); +%% Specify initial parameters +T_sig_len = params.signal.T_sig_len; +T_sig = params.signal.T_sig; +initial_phase = 0; +%% Sampling Frequency +Fs = params.sampling_freq; +%% Number of samples = num*Fs*T_sig +num = params.signal.num; +N = floor(num*T_sig_len*Fs); +timeVec = (0:N-1)*(1/Fs); +%% Min and Max Frequencies Hz +fmin = params.freq(1); +fmax = params.freq(2); +%% Positive Frequency Vector +datalen = N/Fs; +fpos = (0:floor(N/2))*(1/datalen); +%% Initial Time of Arrival and Coalescence Phase +ta = params.ta; +phase = params.phase; +%% Signal to noise ratio of the true signal +snr = params.snr; +%% Constants +c = 3*10^8; +Msolar = 1.989*10^30; +G = 6.6743*10^-11; +%% Mass parameters of the true signal +m1 = params.masses(1); +m2 = params.masses(2); +%% Tau coeffs as phase parameters +if pso.type == "tau" + m1 = m1*Msolar; + m2 = m2*Msolar; + M = m1 + m2; + u = m1*m2/(m1 + m2); + n = u/M; + tau0 = (5/(256*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-5/3))*(1/n); + tau1p5 = (1/8)*(1/fmin)*((G*M*pi*fmin/c^3)^(-2/3))*(1/n); + type = 1; + disp("Tau Space PSO"); +else + type = 0; + disp("Mass Space PSO"); +end +%% Search range of phase coefficients +if type + rmin = [params.rmin_tau(1), params.rmin_tau(2)]; + rmax = [params.rmax_tau(1), params.rmax_tau(2)]; +else + rmin = [params.rmin(1), params.rmin(2)]; + rmax = [params.rmax(1), params.rmax(2)]; +end +%% Number of independent PSO runs +nRuns = pso.nruns; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Optional: Generate Custom Colored Noise +%Generate colored noise using the LIGO Design Sensitivity values using a +%random normal noise realization. Returns, colored noise and two-sided PSD +%vector for postive DFT frequencies, if used Tukey-windowing data is +%recommended to account for any start-up transients. + +% [noise,PSD] = LIGOnoise(N,Fs, params.signal.noise, files.noisefile); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% Pre-processing +[A,avec, phaseDiff] = preprocessing(fmin,fmax, fpos, datalen, N); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% SHAPES and WELCH PSD Comparison Case +%% If used, uncomment lines using PSD and psdVec4norm and comment lines using TF and TFTotal +% %% Generate Final Signal +% S = load(files.datafile); +% dataY = S.dataY; +% %% Load PSD +% E = load(files.psdfile); +% if params.original +% PSD = E.WELCHPSD; +% else +% PSD = E.SHAPESPSD; +% end +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% LIGO HDF5 Datafile Case +%To read conditioned or unconditioned hdf5 data files acquired from GWOSC or GravitySpy. + +dataY = h5read(files.datafile,'/strain/Strain')'; + +%The conditioning transfer function is stored in /strain/condTF by +%convention. If used, comment lines corresponding to vectors PSD and +%psdVec4norm in the following code. +%% If hdf5 file is not conditioned, comment following line +TF = h5read(files.datafile,'/strain/condTF')'; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + + +%% Create entire PSD vector +negFStrt = 1-mod(N,2); +kNyq = floor(N/2)+1; + + +%% Uncomment following lines in the case of SHAPES and WELCH PSD, leave commented otherwise +% psdVec4Norm = [PSD,PSD((kNyq-negFStrt):-1:2)]; +% AbysqrtPSD = A./sqrt(psdVec4Norm); + +%% Uncomment following lines in the case of Conditioned LIGO HDF5 file, leave commented otherwise +TFtotal = [TF, TF((kNyq-negFStrt):-1:2)]; +AbysqrtPSD = A.*TFtotal; + +%% Create General Normalization Factor +% Scalar factor of 1/N is due to Parseval's theorem +dataLen = N; +innProd = (1/dataLen)*(AbysqrtPSD)*AbysqrtPSD'; +genNormfacSqr = real(innProd); +genNormfac = 1/sqrt(genNormfacSqr); + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Optional: Generate and Inject custom CBC signal +%This block of code reads input signal parameters from signal.json and +%creates a custom CBC signal in the time domain that can be injected into +%the data realization. + +% if type +% wavephase = gen2PNwaveform_tau(fpos, ta, phase, fmin, fmax,tau0,tau1p5,datalen, initial_phase, snr, N, avec, genNormfac); +% else +% wavephase = gen2PNwaveform(fpos, ta, phase, fmin, fmax, m1,m2,datalen, initial_phase, snr, N, avec, genNormfac); +% end +% +% wavefourier = A.*wavephase; +% %% Whitening the injected CBC signal to be consistent with strain data +% %% Uncomment following line in the case of Conditioned LIGO HDF5 file, leave commented otherwise +% wavefourier = wavefourier.*TFtotal; + +% %% Uncomment following lines in the case of SHAPES and WELCH PSD, leave commented otherwise +% wavefourier = wavefourier.*psdVec4Norm; + +% wave = ifft(wavefourier); +% +% %% Inject CBC signal into strain data +% dataY = dataY + wave; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% Optional: Tukey-Windowing the data. +%% This to prevent PSO detecting starup transients as signals in the case +%% of custom-made colored noise + +% win = window(@tukeywin, length(dataY)); +% dataY = dataY.*win'; +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% Data Products +fftdataY = fft(dataY); +fftdataY = fftdataY.*A;%Pre-multiply frequency magnitude vector A for optimization + +%% Get FFT of data by total PSD + +%% Uncomment following line in the case of SHAPES and WELCH PSD, leave commented otherwise +%fftdataYbyPSD = fftdataY./psdVec4Norm; + +%% Uncomment following line in the case of Conditioned LIGO HDF5 file, leave commented otherwise +fftdataYbyPSD = fftdataY.*TFtotal; + +dataX = timeVec; + +%% Input Parameters: +inParams = struct('dataX', dataX,... + 'fpos', fpos,... + 'dataY', dataY,... + 'fftdataYbyPSD', fftdataYbyPSD,... + 'frange', [fmin,fmax],... + 'datalen',datalen,..., + 'initial_phase', initial_phase,... + 'N', N,... + 'A', A,... + 'phaseDiff', phaseDiff,... + 'normfac', genNormfac,... + 'avec', avec,... + 'T_sig', T_sig,... + 'rmin',rmin,... + 'rmax',rmax,... + 'Fs',Fs); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +maxSteps = pso.maxSteps; +if type + original_fitVal = -1*mfgw_tau([tau0, tau1p5], inParams); + outStruct = crcbgwpso_tau(inParams,struct('maxSteps',maxSteps),nRuns,Fs); + bestFitVal = -1*outStruct.bestFitness; +else + original_fitVal = -1*mfgw_mass([m1, m2], inParams); + outStruct = crcbgwpso_mass(inParams,struct('maxSteps',maxSteps),nRuns,Fs); + bestFitVal = -1*outStruct.bestFitness; +end + +%% Optional: Uncomment to save output of PSO in a .mat file +% save(files.output_struct_location,'outStruct'); + +%% Plots, uncomment all saveas() commands to save figures in custom directories +% psooutput(outStruct,maxSteps,dataX, nRuns, dataY, tau0, tau1p5, m1, m2, snr, ta, phase, original_fitVal, bestFitVal, c, G, type); +figure; +hold on; +plot(dataX,dataY,'.'); +% plot(dataX,wave,'r'); +% for lpruns = 1:nRuns +% plot(dataX,outStruct.allRunsOutput(lpruns).estSig,'Color',[51,255,153]/255,'LineWidth',4.0); +% end +plot(dataX,outStruct.bestSig,'Color',[76,153,0]/255,'LineWidth',2.0); +% legend('Data','Signal',... +% ['Estimated signal: ',num2str(nRuns),' runs'],... +% 'Estimated signal: Best run'); +legend('Data','Signal',... + 'Estimated signal: Best run'); +% saveas(gcf,files.psoresultplot); +hold off; + +figure; +iterVec = linspace(1,maxSteps,maxSteps); +hold on; +for lpruns = 1:nRuns + plot(iterVec,outStruct.allRunsOutput(lpruns).allBestFit, 'DisplayName',num2str(lpruns)); +end +title("Best Fitness Values for All Runs"); +xlabel("Iteration"); +ylabel("Best Fitness Value"); +legend; +% saveas(gcf,files.bestfitplot); +hold off; + +if type + figure; + hold on; + for lpruns = 1:nRuns + rVec = s2rv(outStruct.allRunsOutput(lpruns).allBestLoc,inParams); + plot(rVec(:,1),rVec(:,2),'DisplayName',num2str(lpruns)); + end + scatter(tau0,tau1p5,140,'red','filled','D','DisplayName','Original Parameters'); + title("Best Parameter Values for All Runs"); + xlabel("\tau_0"); + ylabel("\tau_{1.5}"); + legend; + boundary_plot; +% saveas(gcf,files.bestlocplot); + hold off; + + + + t0 = outStruct.bestQcCoefs(1); + t1p5 = outStruct.bestQcCoefs(2); + est_M = (5/(32*fmin))*(t1p5/(pi*pi*t0))*(c^3/G); + est_u = (1/(16*fmin*fmin))*(5/(4*pi^4*t0*t1p5^2))^(1/3)*(c^3/G); + + est_m1 = (est_M - sqrt(est_M^2 - 4*est_u*est_M))/2; + est_m2 = (est_M + sqrt(est_M^2 - 4*est_u*est_M))/2; + + %% This will display parameters given through signal.json and PSO-estimated parameters + %% Uncomment Original parameter display command if needed + disp(['Original parameters: tau0= ',num2str(tau0),... + '; tau1p5= ',num2str(tau1p5),... + '; m1= ', num2str(m1/Msolar),... + '; m2= ', num2str(m2/Msolar),... + '; A = ',num2str(snr),... + '; phi = ',num2str(phase),... + '; t_a = ',num2str(ta),... + '; FitVal = ',num2str(original_fitVal)]); + + disp(['Estimated parameters: tau0=',num2str(outStruct.bestQcCoefs(1)),... + '; tau1p5=',num2str(outStruct.bestQcCoefs(2)),... + '; m1= ', num2str(est_m1/Msolar),... + '; m2= ', num2str(est_m2/Msolar),... + '; A = ',num2str(outStruct.bestAmp),... + '; phi = ',num2str(outStruct.bestPhase),... + '; t_a = ',num2str(outStruct.bestTime),... + '; FitVal = ',num2str(bestFitVal)]); +else + figure; + hold on; + for lpruns = 1:nRuns + rVec = s2rv(outStruct.allRunsOutput(lpruns).allBestLoc,inParams); + plot(rVec(:,1),rVec(:,2),'DisplayName',num2str(lpruns)); + end + scatter(m1,m2,140,'red','filled','D','DisplayName','Original Parameters'); + title("Best Parameter Values for All Runs"); + xlabel("m_1"); + ylabel("m_2"); + legend; +% saveas(gcf,files.bestlocplot); + hold off; + + %% This will display parameters given through signal.json and PSO-estimated parameters + %% Uncomment Original parameter display command if needed + disp(['Original parameters: m1= ',num2str(m1),... + '; m2= ',num2str(m2),... + '; A = ',num2str(snr),... + '; phi = ',num2str(phase),... + '; t_a = ',num2str(ta),... + '; FitVal = ',num2str(original_fitVal)]); + + disp(['Estimated parameters: m1=',num2str(outStruct.bestQcCoefs(1)),... + '; m2=',num2str(outStruct.bestQcCoefs(2)),... + '; A = ',num2str(outStruct.bestAmp),... + '; phi = ',num2str(outStruct.bestPhase),... + '; t_a = ',num2str(outStruct.bestTime),... + '; FitVal = ',num2str(bestFitVal)]); +end +end diff --git a/Misc/NormalizationCheck/segdatacond.m b/Misc/NormalizationCheck/segdatacond.m new file mode 100644 index 0000000..740e926 --- /dev/null +++ b/Misc/NormalizationCheck/segdatacond.m @@ -0,0 +1,48 @@ +function [whtndfiltdata, whtndstd, TFtotal]=segdatacond(segdata, PSD, sampFreq, tidxs) + +fmin = 30; +% sampfreq = 4096; + +filtdata = segdata; + +datalen = length(filtdata)/sampFreq; + +N = length(filtdata); +Tsig = N/sampFreq; +timeVec = (0:N-1)*(1/sampFreq); + +%Highpass Filter the data and handle NaN-regions +% [filtdata, PSD, nanchunk_start_idxs, nanchunk_end_idxs] = nanhandle(data, sampFreq, fmin, tstart, seglen, winlen); + +%Create Transfer Function +TF = 1./sqrt(PSD); + +%Set Values of TF below fmin to 0 +TF(1:fmin*datalen) = 0; + +%Create entire Transfer Function vector +negFStrt = 1-mod(N,2); +kNyq = floor(N/2)+1; + +TFtotal = [TF, TF((kNyq-negFStrt):-1:2)]; + +%Whiten Filter Data with Transfer Function +%Window data before FFT with a Tukey-window with a 0.5 seconds rolloff + +rolloff = 0.5; %Roll-off in seconds +winfiltdata = filtdata.*tukeywin(length(filtdata), rolloff*sampFreq/N)'; + +fftfiltdata = fft(winfiltdata); + +whtndfftfiltdata = fftfiltdata.*TFtotal; + +whtndfiltdata = ifft(whtndfftfiltdata); + + +%Divide by variance of whitened strain so that final whitened vector has +%unit strain +tstart = tidxs(1); +tend = tidxs(2); +whtndstd = std(whtndfiltdata(tstart: tend)); +whtndfiltdata = whtndfiltdata/sqrt(sampFreq); +% whtndfiltdata = whtndfiltdata/whtndstd; diff --git a/Misc/NormalizationCheck/signal.json b/Misc/NormalizationCheck/signal.json new file mode 100644 index 0000000..4e7b756 --- /dev/null +++ b/Misc/NormalizationCheck/signal.json @@ -0,0 +1,27 @@ +{ +"sampling_freq": 4096, +"ta": 900, +"snr": 15, +"phase": 0, +"masses": [2, 2], +"rmin": [1, 1], +"rmax": [30, 30], +"rmin_tau": [0, 0], +"rmax_tau": [90, 2], +"freq": [30, 700], +"original": 0, +"fitnessimage": { +"startidx_tau0": 11, +"startidx_tau1p5": 11, +"ngridpoints": 2000, +"nsectors": 4, +"tau0_range": [0.1, 70], +"tau1p5_range": [0.1, 2] +}, +"signal": { +"T_sig_len": 2048, +"T_sig": 54, +"num": 1, +"noise": 1 +} +} diff --git a/Misc/NormalizationCheck/signal_injection_codes.zip b/Misc/NormalizationCheck/signal_injection_codes.zip new file mode 100644 index 0000000..d6b67d7 Binary files /dev/null and b/Misc/NormalizationCheck/signal_injection_codes.zip differ diff --git a/Misc/NormalizationCheck/statgaussnoisegen.m b/Misc/NormalizationCheck/statgaussnoisegen.m new file mode 100644 index 0000000..4816590 --- /dev/null +++ b/Misc/NormalizationCheck/statgaussnoisegen.m @@ -0,0 +1,46 @@ +function outNoise = statgaussnoisegen(nSamples,psdVals,fltrOrdr,sampFreq, noise_num, noisefile) +%Generate a realization of stationary Gaussian noise with given 2-sided PSD +%Y = STATGAUSSNOISEGEN(N,PSD,fltrOrdr,Fs, noise_num, noisefile) +%Generates a realization Y of stationary gaussian noise with a target +%2-sided power spectral density given by PSD. +% Fs is the sampling frequency +%of Y. +% PSD is an M-by-2 matrix containing frequencies and the corresponding +%PSD values in the first and second columns respectively. +% The frequencies +%must start from 0 and end at Fs/2. +% The order of the FIR filter to be usedis given by fltrOrdr. +%(Optional) noise_num = noise realization number from a pre-created noise realizations file +%(Optional) noisefile = pre-created noise realizations filename + +%Soumya D. Mohanty, Mar 2019 + +% Added Tukey-windowing option and option to load noise realization from +% custom pre-made noise files + +% Raghav Girgaonkar, Apr 2023 + +%% Design FIR filter with T(f)= square root of target PSD +freqVec = psdVals(:,1); +sqrtPSD = sqrt(psdVals(:,2)); +b = fir2(fltrOrdr,freqVec/(sampFreq/2),sqrtPSD,tukeywin(fltrOrdr+1)); + +%% Generate a WGN realization and pass it through the designed filter + +%% Uncomment following block of code if using pre-created normal noise file +%% leave commented otherwise + +% noise = load(noisefile); +% inNoise_t = noise.wgn(noise_num,1:end); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +%% Uncomment following block of code if creating noise vector +%% leave commented otherwise +% (Comment out the line below if new realizations of WGN are needed in each run of this script) +% rng('default'); + +inNoise_t = randn(1,nSamples); +inNoise = [randn(1,fltrOrdr), inNoise_t, randn(1,fltrOrdr)]; +outNoise = sqrt(sampFreq)*fftfilt(b,inNoise); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + diff --git a/Misc/NormalizationCheck/test_matchedfiltering.m b/Misc/NormalizationCheck/test_matchedfiltering.m new file mode 100644 index 0000000..c527a25 --- /dev/null +++ b/Misc/NormalizationCheck/test_matchedfiltering.m @@ -0,0 +1,314 @@ +% Define Parameters +sampFreq = 4096; +datalen = 512; +N = datalen*sampFreq; +frange = [30,700]; +fpos = (0:floor(N/2))*(1/datalen); +timeVec = (0:(datalen*sampFreq - 1))/sampFreq; + +%Create Colored Gaussian Noise +[noise, PSD] = LIGOnoise(N, sampFreq, 1, 'sample'); + +%Create Total PSD +negFStrt = 1-mod(N,2); +kNyq = floor(N/2)+1; + +PSDtotal = [PSD, PSD((kNyq-negFStrt):-1:2)]; + +%Generate signal normalized to specified SNR + chirptimes = [29.6373, 1.1045]; + ta = 138; + initial_phase = 0; + phase = 0; + snr = 30; + + [A,avec, phaseDiff] = preprocessing(frange(1),frange(2),fpos, datalen, N); + + signal = gen2PNtemplate_tau(fpos, ta, phase, frange(1), frange(2),chirptimes(1), chirptimes(2),datalen,initial_phase,snr, N, A, avec, PSDtotal); + + signal = signal*sqrt(sampFreq); +% signal = signal/norm(signal); + + figure; + plot(timeVec,signal); + + %Test 1: Inject signal directly into strain data + +straindata = noise + signal; + + % Input Parameters Structure: +inParams = struct('fpos', fpos,... + 'dataY', straindata,... + 'frange', frange,... + 'datalen',datalen,..., + 'initial_phase', initial_phase,... + 'N', N,... + 'A', A,... + 'phaseDiff', phaseDiff,... + 'avec', avec,... + 'T_sig', 54,... + 'PSDtotal',PSDtotal,... + 'Fs',sampFreq); + +[mf1,mf2,max_val,max_arg] = mfgw_tau_new(chirptimes, inParams); +mftimeseries = sqrt(mf1(1:end - 54*sampFreq).^2 + mf2(1:end - 54*sampFreq).^2); + +figure; +plot(timeVec,mf1); title('MF1 Test1'); xlabel('t'); + +figure; +plot(timeVec,mf2); title('MF2 Test1'); xlabel('t'); + +figure; +plot(timeVec(1:end - 54*sampFreq),mftimeseries); title('MFTimeseries Test1'); xlabel('t'); + + +%Test 2: Inject signal into strain data and whiten together without +%normalizing whitened strain to unit variance + +[whtndseg,wstdv, TFtotal] = segdatacond(straindata, PSD, sampFreq, [1,20*sampFreq]); + +% whtndseg = whtndseg*wstdv; + +% Input Parameters Structure: +inParams = struct('fpos', fpos,... + 'dataY', whtndseg,... + 'frange', frange,... + 'datalen',datalen,..., + 'initial_phase', initial_phase,... + 'N', N,... + 'A', A,... + 'phaseDiff', phaseDiff,... + 'avec', avec,... + 'T_sig', 54,... + 'PSDtotal',PSDtotal,... + 'TFtotal',TFtotal,... + 'Fs',sampFreq); + +[mf1,mf2,max_val,max_arg] = mfgw_tau_whtnd_og(chirptimes, inParams); +mftimeseries = sqrt(mf1(1:end - 54*sampFreq).^2 + mf2(1:end - 54*sampFreq).^2); + +figure; +plot(timeVec,mf1); title('MF1 Test2'); xlabel('t'); + +figure; +plot(timeVec,mf2); title('MF2 Test2'); xlabel('t'); + +figure; +plot(timeVec(1:end - 54*sampFreq),mftimeseries); title('MFTimeseries Test2'); xlabel('t'); + + +%Test 3: Inject signal into strain data and whiten together with +%normalizing whitened strain to unit variance + +[whtndseg,~, TFtotal] = segdatacond(straindata, PSD, sampFreq, [1,20*sampFreq]); + + +% Input Parameters Structure: +inParams = struct('fpos', fpos,... + 'dataY', whtndseg,... + 'frange', frange,... + 'datalen',datalen,..., + 'initial_phase', initial_phase,... + 'N', N,... + 'A', A,... + 'phaseDiff', phaseDiff,... + 'avec', avec,... + 'T_sig', 54,... + 'PSDtotal',PSDtotal,... + 'TFtotal',TFtotal,... + 'Fs',sampFreq); + +[mf1,mf2,max_val,max_arg] = mfgw_tau_whtnd_og(chirptimes, inParams); +mftimeseries = sqrt(mf1(1:end - 54*sampFreq).^2 + mf2(1:end - 54*sampFreq).^2); + +figure; +plot(timeVec,mf1); title('MF1 Test3'); xlabel('t'); + +figure; +plot(timeVec,mf2); title('MF2 Test3'); xlabel('t'); + +figure; +plot(timeVec(1:end - 54*sampFreq),mftimeseries); title('MFTimeseries Test3'); xlabel('t'); + + +%Test 4: Inject a whitened signal into a whitened data segment which has +%been normalized to unit variance. + +%First generate normalized whiten data realization +[whtndseg,~, TFtotal] = segdatacond(noise, PSD, sampFreq, [1,20*sampFreq]); + + whtndseg = whtndseg; + +%Generate whitened signal +whtndsig = ifft(fft(signal).*(TFtotal)); + +%Add whitened signal to normalized whitened data +whtndseg = whtndseg + whtndsig; + +% Input Parameters Structure: +inParams = struct('fpos', fpos,... + 'dataY', whtndseg,... + 'frange', frange,... + 'datalen',datalen,..., + 'initial_phase', initial_phase,... + 'N', N,... + 'A', A,... + 'phaseDiff', phaseDiff,... + 'avec', avec,... + 'T_sig', 54,... + 'PSDtotal',PSDtotal,... + 'TFtotal',TFtotal,... + 'Fs',sampFreq); + +[mf1,mf2,max_val,max_arg] = mfgw_tau_whtnd_og(chirptimes, inParams); +mftimeseries = sqrt(mf1(1:end - 54*sampFreq).^2 + mf2(1:end - 54*sampFreq).^2); + +figure; +plot(timeVec,mf1); title('MF1 Test4'); xlabel('t'); + +figure; +plot(timeVec,mf2); title('MF2 Test4'); xlabel('t'); + +figure; +plot(timeVec(1:end - 54*sampFreq),mftimeseries); title('MFTimeseries Test4'); xlabel('t'); + +%Test 5, add injection directly into strain data, but scaled by the stdv of +%the whitened data. + +%Get stdv +[~,wstdv,~] = segdatacond(noise, PSD, sampFreq, [1,200*sampFreq]); + +% Create signal and scale + +sigscaled = signal*wstdv; + +%Add to strain + +straindata = noise + sigscaled; + + % Input Parameters Structure: +inParams = struct('fpos', fpos,... + 'dataY', straindata,... + 'frange', frange,... + 'datalen',datalen,..., + 'initial_phase', initial_phase,... + 'N', N,... + 'A', A,... + 'phaseDiff', phaseDiff,... + 'avec', avec,... + 'T_sig', 54,... + 'PSDtotal',PSDtotal,... + 'Fs',sampFreq); + +[mf1,mf2,max_val,max_arg] = mfgw_tau_new(chirptimes, inParams); +mftimeseries = sqrt(mf1(1:end - 54*sampFreq).^2 + mf2(1:end - 54*sampFreq).^2); + +figure; +plot(timeVec,mf1); title('MF1 Test5'); xlabel('t'); + +figure; +plot(timeVec,mf2); title('MF2 Test5'); xlabel('t'); + +figure; +plot(timeVec(1:end - 54*sampFreq),mftimeseries); title('MFTimeseries Test5'); xlabel('t'); + + + +%Test 6, matched filtering on whitened data normalized to unit variance: + + +%First generate normalized whiten data realization +[whtndseg,~, TFtotal] = segdatacond(noise, PSD, sampFreq, [1,datalen*sampFreq]); + +% Input Parameters Structure: +inParams = struct('fpos', fpos,... + 'dataY', whtndseg,... + 'frange', frange,... + 'datalen',datalen,..., + 'initial_phase', initial_phase,... + 'N', N,... + 'A', A,... + 'phaseDiff', phaseDiff,... + 'avec', avec,... + 'T_sig', 54,... + 'PSDtotal',PSDtotal,... + 'TFtotal',TFtotal,... + 'Fs',sampFreq); + +[mf1,mf2,max_val,max_arg] = mfgw_tau_whtnd_og(chirptimes, inParams); +mftimeseries = sqrt(mf1(1:end - 54*sampFreq).^2 + mf2(1:end - 54*sampFreq).^2); + +figure; +plot(timeVec,mf1); title('MF1 Test6'); xlabel('t'); + +figure; +plot(timeVec,mf2); title('MF2 Test6'); xlabel('t'); + +figure; +plot(timeVec(1:end - 54*sampFreq),mftimeseries); title('MFTimeseries Test6'); xlabel('t'); + +%Test 7: Inject signal directly into the strain data scaled by the whtnd +%stdv + +[~,wstdv, ~] = segdatacond(noise, PSD, sampFreq, [1,datalen*sampFreq]); + +straindata = noise + wstdv*signal; + +[whtndseg,~, TFtotal] = segdatacond(straindata, PSD, sampFreq, [1,20*sampFreq]); + + % Input Parameters Structure: +inParams = struct('fpos', fpos,... + 'dataY', whtndseg,... + 'frange', frange,... + 'datalen',datalen,..., + 'initial_phase', initial_phase,... + 'N', N,... + 'A', A,... + 'phaseDiff', phaseDiff,... + 'avec', avec,... + 'T_sig', 54,... + 'PSDtotal',PSDtotal,... + 'TFtotal',TFtotal,... + 'Fs',sampFreq); + +[mf1,mf2,max_val,max_arg] = mfgw_tau_whtnd_og(chirptimes, inParams); +mftimeseries = sqrt(mf1(1:end - 54*sampFreq).^2 + mf2(1:end - 54*sampFreq).^2); + +figure; +plot(timeVec,mf1); title('MF1 Test7'); xlabel('t'); + +figure; +plot(timeVec,mf2); title('MF2 Test7'); xlabel('t'); + +figure; +plot(timeVec(1:end - 54*sampFreq),mftimeseries); title('MFTimeseries Test7'); xlabel('t'); + + +%Test 8: Repeat Test 7 100 times and get histogram of recovered SNR + +snrs = []; +for i = 1:100 + disp(i); + [noise, ~] = LIGOnoise(N, sampFreq, 1, 'sample'); + [~,wstdv, ~] = segdatacond(noise, PSD, sampFreq, [1,datalen*sampFreq]); + + straindata = noise + wstdv*signal; + + [whtndseg,~, TFtotal] = segdatacond(straindata, PSD, sampFreq, [1,20*sampFreq]); + inParams.dataY = whtndseg; + [~,~,max_val,~] = mfgw_tau_whtnd_og(chirptimes, inParams); + snrs = [snrs, sqrt(-max_val)]; +end + +figure; +histogram(snrs, 50); +hold on; +xline(snr); +hold off; + + + + + + diff --git a/Misc/NormalizationCheck/waveform.m b/Misc/NormalizationCheck/waveform.m new file mode 100644 index 0000000..623eb8a --- /dev/null +++ b/Misc/NormalizationCheck/waveform.m @@ -0,0 +1,73 @@ +function fwave = waveform(fvec, t, phase, fmin, fmax,m1, m2,datalen,initial_phase, avec) +%Function to create Restricted 2PN Waveform in Fourier Domain +%Input: fvec = Row vector of positive DTFT frequencies +% t = Time of arrival in seconds +% phase = Phase of Coalescence +% fmin, fmax = minimum and maximum frequency cutoffs for the waveform +% m1, m2 = Component masses +% datalen = Length of data realization in seconds +% avec = Pre-calculated alpha term vector +%Output: fwave = Phase term of waveform in the Fourier domain + +%Raghav Girgaonkar, April 2023 + +%% Constants +c = 3*10^8; +Msolar = 1.989*10^30; +G = 6.6743*10^-11; +%% Calculate Mass Terms +m1 = m1*Msolar; +m2 = m2*Msolar; +M = (m1 + m2); +u = m1*m2/M; +n = u/M; +%% Calculate Chirp Times +tau0 = (5/(256*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-5/3))*(1/n); + +tau1 = (5/(192*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-1))*(1/n)*((743/336)+ (11*n/4)); + +tau1p5 = (1/8)*(1/fmin)*((G*M*pi*fmin/c^3)^(-2/3))*(1/n); + +tau2 = (5/(128*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-1/3))*(1/n)*((3058673/1016064) + (5429*n/1008) + (617*n*n/144)); + +%% Alpha Terms +alpha0 = 2*pi*fmin*(3*tau0/5); + +alpha1 = 0; + +alpha2 = 2*pi*fmin*tau1; + +alpha3 = -2*pi*fmin*(3*tau1p5/2); + +alpha4 = 2*pi*fmin*3*tau2; + +alphaTerm = zeros(size(fvec)); + +alphaTerm(2:end) = alpha0*avec(1,:)... ++ alpha1*avec(2,:)... ++ alpha2*avec(3,:)... ++ alpha3*avec(4,:)... ++ alpha4*avec(5,:); + +F = 2*pi*fvec*(tau0 + tau1 - tau1p5 + tau2); + +%% Final Phase Term + +Psi = 2*pi*t*fvec - phase - pi/4 + alphaTerm + F + initial_phase; + + +%% Final Expression +P = exp(-1j*Psi); +fwave = P; + +%% Cut between fmin and fmax +min_index = floor(datalen*fmin) + 1; +max_index = floor(datalen*fmax) + 1; + +fwave(1:min_index-1) = 0; +fwave(max_index + 1: end) = 0; + + + + + diff --git a/Misc/NormalizationCheck/waveform_tau.m b/Misc/NormalizationCheck/waveform_tau.m new file mode 100644 index 0000000..6003f6e --- /dev/null +++ b/Misc/NormalizationCheck/waveform_tau.m @@ -0,0 +1,80 @@ +function fwave = waveform_tau(fvec, t, phase, fmin, fmax,tau0, tau1p5,datalen,initial_phase,avec) +%Function to create Restricted 2PN Waveform in Fourier Domain +%Input: fvec = Row vector of positive DTFT frequencies +% t = Time of arrival in seconds +% phase = Phase of Coalescence +% fmin, fmax = minimum and maximum frequency cutoffs for the waveform +% tau0, tau1.5 = Chirp-time parameters +% datalen = Length of data realization in seconds +% avec = Pre-calculated alpha term vector +%Output: fwave = Phase term of waveform in the Fourier domain + +%Raghav Girgaonkar, April 2023 + +%% Constants +% c = 299792458; +c = 3*10^8; +% Msolar = 1.989*10^30; +G = 6.6743*10^-11; +%% Calculate Mass Terms from Tau0 and Tau1.5 +M = (5/(32*fmin))*(tau1p5/(pi*pi*tau0))*(c^3/G); +if tau0 < 0 && tau1p5 < 0 + u = -(1/(16*fmin*fmin))*(5/(4*pi^4*abs(tau0)*abs(tau1p5)^2))^(1/3)*(c^3/G); +else + u = (1/(16*fmin*fmin))*(5/(4*pi^4*tau0*tau1p5^2))^(1/3)*(c^3/G); +end +% m1 = m1*Msolar; +% m2 = m2*Msolar; + +n = u/M; +%% Calculate Chirp Times + +% tau0 = (5/(256*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-5/3))*(1/n); + +tau1 = (5/(192*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-1))*(1/n)*((743/336)+ (11*n/4)); + +% tau1p5 = (1/8)*(1/fmin)*((G*M*pi*fmin/c^3)^(-2/3))*(1/n); + +tau2 = (5/(128*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-1/3))*(1/n)*((3058673/1016064) + (5429*n/1008) + (617*n*n/144)); + +%Alpha Terms +alpha0 = 2*pi*fmin*(3*tau0/5); + +alpha1 = 0; + +alpha2 = 2*pi*fmin*tau1; + +alpha3 = -2*pi*fmin*(3*tau1p5/2); + +alpha4 = 2*pi*fmin*3*tau2; + +alphaTerm = zeros(size(fvec)); + +alphaTerm(2:end) = alpha0*avec(1,:)... ++ alpha1*avec(2,:)... ++ alpha2*avec(3,:)... ++ alpha3*avec(4,:)... ++ alpha4*avec(5,:); + +F = 2*pi*fvec*(tau0 + tau1 - tau1p5 + tau2); + +%% Final Phase Term + +Psi = 2*pi*t*fvec - phase - pi/4 + alphaTerm + F + initial_phase; + + +%% Final Expression +P = exp(-1j*Psi); +fwave = P; + +%% Cut between fmin and fmax +min_index = floor(datalen*fmin) + 1; +max_index = floor(datalen*fmax) + 1; + +fwave(1:min_index-1) = 0; +fwave(max_index + 1: end) = 0; + + + + + diff --git a/Misc/NormalizationCheck/waveform_tau_negative.m b/Misc/NormalizationCheck/waveform_tau_negative.m new file mode 100644 index 0000000..9793846 --- /dev/null +++ b/Misc/NormalizationCheck/waveform_tau_negative.m @@ -0,0 +1,102 @@ +function fwave = waveform_tau_negative(fvec, t, phase, fmin, fmax,tau0, tau1p5,datalen,initial_phase,avec) +%Function to create Restricted 2PN Waveform in Fourier Domain in the +%Negative Chirp Time space +%Input: fvec = Row vector of positive DTFT frequencies +% t = Time of arrival in seconds +% phase = Phase of Coalescence +% fmin, fmax = minimum and maximum frequency cutoffs for the waveform +% tau0, tau1.5 = Chirp-time parameters +% datalen = Length of data realization in seconds +% avec = Pre-calculated alpha term vector +%Output: fwave = Phase term of waveform in the Fourier domain + +%Raghav Girgaonkar, April 2023 + +%% Constants +c = 3*10^8; +% Msolar = 1.989*10^30; +G = 6.6743*10^-11; + +%% Check for positive chirp-length and swap tau0 and tau1.5 values if yes + + M = (5/(32*fmin))*(tau1p5/(pi*pi*tau0))*(c^3/G); + + if tau0 < 0 && tau1p5 < 0 + u = 1/(16*fmin*fmin)*(5/(4*pi^4*abs(tau0)*abs(tau1p5)^2))^(1/3)*(c^3/G); + else + u = (1/(16*fmin*fmin))*(5/(4*pi^4*tau0*tau1p5^2))^(1/3)*(c^3/G); + end + n = u/M; + %% Calculate Chirp Times + tau1 = -(5/(192*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-1))*(1/n)*((743/336)+ (11*n/4)); + tau2 = -(5/(128*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-1/3))*(1/n)*((3058673/1016064) + (5429*n/1008) + (617*n*n/144)); + + if tau0 - tau1p5 + tau1 + tau2 > 0 + temp = tau0; + tau0 = tau1p5; + tau1p5 = temp; + end + + +%% Calculate Mass Terms from Tau0 and Tau1.5 +M = (5/(32*fmin))*(tau1p5/(pi*pi*tau0))*(c^3/G); +if tau0 < 0 && tau1p5 < 0 + u = 1/(16*fmin*fmin)*(5/(4*pi^4*abs(tau0)*abs(tau1p5)^2))^(1/3)*(c^3/G); +else + u = (1/(16*fmin*fmin))*(5/(4*pi^4*tau0*tau1p5^2))^(1/3)*(c^3/G); +end +% m1 = m1*Msolar; +% m2 = m2*Msolar; + +n = u/M; +%% Calculate Chirp Times + +% tau0 = (5/(256*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-5/3))*(1/n); + +tau1 = -(5/(192*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-1))*(1/n)*((743/336)+ (11*n/4)); + +% tau1p5 = (1/8)*(1/fmin)*((G*M*pi*fmin/c^3)^(-2/3))*(1/n); + +tau2 = -(5/(128*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-1/3))*(1/n)*((3058673/1016064) + (5429*n/1008) + (617*n*n/144)); + +%Alpha Terms +alpha0 = 2*pi*fmin*(3*tau0/5); + +alpha1 = 0; + +alpha2 = 2*pi*fmin*tau1; + +alpha3 = -2*pi*fmin*(3*tau1p5/2); + +alpha4 = 2*pi*fmin*3*tau2; + +alphaTerm = zeros(size(fvec)); + +alphaTerm(2:end) = alpha0*avec(1,:)... ++ alpha1*avec(2,:)... ++ alpha2*avec(3,:)... ++ alpha3*avec(4,:)... ++ alpha4*avec(5,:); + +% F = 0*2*pi*fvec*(abs(tau0) + abs(tau1) - abs(tau1p5) + abs(tau2)); + +%% Final Phase Term + +Psi = 2*pi*t*fvec - phase - pi/4 + alphaTerm + initial_phase; + + +%% Final Expression +P = exp(-1j*Psi); +fwave = P; + +%% Cut between fmin and fmax +min_index = floor(datalen*fmin) + 1; +max_index = floor(datalen*fmax) + 1; + +fwave(1:min_index-1) = 0; +fwave(max_index + 1: end) = 0; + + + + + diff --git a/Note_On_Normalization.pdf b/Note_On_Normalization.pdf new file mode 100644 index 0000000..f72e691 Binary files /dev/null and b/Note_On_Normalization.pdf differ diff --git a/OptimizedPSOCodes/launcherscript.m b/OptimizedPSOCodes/launcherscript.m new file mode 100755 index 0000000..046535d --- /dev/null +++ b/OptimizedPSOCodes/launcherscript.m @@ -0,0 +1,157 @@ +function []=launcherscript(launcherparamfilename) +%Script to run PSO-based Matched Filtering on a data segment +%Inputs: segment number, txt file with filenames, txt file with segment +%details + +%Load Files +fname = launcherparamfilename; +str = fileread(fname); +launcherparams = jsondecode(str); +hdf5filenames = launcherparams.hdf5filenames; +segfile = launcherparams.segfile; +psdfile = launcherparams.outfilename; +sampFreq = launcherparams.sampFreq; +segnum = launcherparams.segnum; +threshold = launcherparams.threshold; + + +%Load Data vector from file(s) +%Get Segment Data +segments = textread(segfile, '%s', 'delimiter', '\n'); +Seg = split(segments{segnum}); +segstart = str2num(Seg{2}); +segend = str2num(Seg{3}); +filenums = str2num(Seg{4}); + +%Get File data and load segment +if length(filenums) == 1 + files = textread(hdf5filenames, '%s', 'delimiter', '\n'); + fileStr = split(files{filenums}); + filename = fileStr{1}; + GPSstart = str2num(fileStr{2}); + + %Get segment indices + startidx = sampFreq*(segstart- GPSstart); + endidx = sampFreq*(segend- GPSstart); + + %Load data from file and create segment data vector + filedata = h5read(filename, '/strain/Strain')'; + segment = filedata(startidx+1:endidx); + psdnum = filenums; +else + segment = []; + psdnum = filenums(1); + for i = 1:length(filenums) + files = textread(hdf5filenames, '%s', 'delimiter', '\n'); + fileStr = split(files{filenums(i)}); + filename = fileStr{1}; + GPSstart = str2num(fileStr{2}); + GPSend = str2num(fileStr{3}); + + if GPSend > segstart && GPSstart < segstart && segend > GPSend + %Get segment indices + startidx = sampFreq*(segstart- GPSstart); + endidx = sampFreq*(GPSend -GPSstart); + %Load data from file and create segment data vector + filedata = h5read(filename, '/strain/Strain')'; + seg = filedata(startidx:endidx); + segment = [segment, seg]; + end + if GPSstart < segend && GPSend > segend && segstart < GPSstart + %Get segment indices + startidx = 1; + endidx = sampFreq*(segend - GPSstart); + %Load data from file and create segment data vector + filedata = h5read(filename, '/strain/Strain')'; + seg = filedata(startidx:endidx-1); + segment = [segment, seg]; + end + end +end + +Seglen = length(segment)/sampFreq; + +%{ +signalinjectionsegs = [7,9,13,16,17,18]; + +%Inject signals +fmin = 30; +fmax = 700; +frange = [fmin,fmax]; +initial_phase = 0; phase = 0; +siglen = Seglen; +if sum(ismember(signalinjectionsegs,segnum)) == 1 + %Draw parameters from uniform distribution + m1 = unifrnd(1.4,15); + m2 = unifrnd(1.4,15); + masses = [m1,m2]; + %Estimate the corresponding chirptime parameters + m1_val = m1*Msolar; + m2_val = m2*Msolar; + M = m1_val + m2_val; + u = m1_val*m2_val/(m1_val + m2_val); + n = u/M; + tau0 = (5/(256*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-5/3))*(1/n); + tau1p5 = (1/8)*(1/fmin)*((G*M*pi*fmin/c^3)^(-2/3))*(1/n); + + r = unifrnd(30,100); + mfac = unifrnd(2,3.2); + ta = unifrnd(100,200); + + %Create signal + signal = createsignal(siglen, frange, sampFreq, masses, r, initial_phase, phase, ta, mfac); + + %Record simulated strain signal in file + sigparams = [segnum, tau0, tau1p5, m1, m2, ta, r, mfac]; + fileID = fopen('injectedsigs.txt','a'); + fprintf(fileID,'%d\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n',sigparams); + fclose(fileID); + %Add signal to segment + segment = segment + signal; +end +%} + +%Load Training Segment PSD +S = load(psdfile); +segPSDs = S.segPSDs; +trainidxs = S.trainidxs; +PSD = segPSDs{segnum}; +trainidxs = trainidxs{segnum}; + + +%Condition Data +rolloff = 0.125; +fmin = 30; +% paddedsegment = [zeros(1,sampFreq*rolloff), segment, zeros(1,sampFreq*rolloff)]; +%Window and Highpass filter +segwin = segment.*tukeywin(length(segment),rolloff*sampFreq/length(segment))'; +seghpass = highpass(segwin, fmin, sampFreq, ImpulseResponse="iir",Steepness=0.95); +filtsegment = seghpass; + + +%Whiten Segment +timeVec = (0:(length(filtsegment) - 1))/sampFreq; +[whtndseg,TFtotal] = segdatacond(filtsegment, PSD, sampFreq); + + +% Add pre-created, whitened custom signals in specific segments +%S = load('GVSsiginjections_final.mat'); sigsegnum = S.sigsegnums; signals = S.signals; +% S = load('GVSsiginjections_highsnr.mat'); sigsegnum = S.siginjsegs; signals = S.signals; +% index = find(sigsegnum == segnum); +% +% if index +% signal = signals(index,:); +% whtndseg = whtndseg + signal; +% end +% +% %Special modification for negative chirp-time space run on segment 878 which contains LIGO event GW190707_093326: Replaces the glitch near 347.8 seconds with WGN +% +% if segnum == 878 +% whtndseg(347*sampFreq:348*sampFreq) = randn(1, length(whtndseg(347*sampFreq:348*sampFreq))); +% end + + +%Run PSO based matched-filtering on data vector +paramfile = 'allparamfiles.json'; +runpso(segment, whtndseg, TFtotal, paramfile,segnum,Seglen); +end diff --git a/OptimizedPSOCodes/runpso.m b/OptimizedPSOCodes/runpso.m new file mode 100755 index 0000000..40a0148 --- /dev/null +++ b/OptimizedPSOCodes/runpso.m @@ -0,0 +1,271 @@ +function [] = runpso(segdata, whtnddata,TFtotal, filename, segnum,Seglen) +%Script to run Chirp-time or Mass Space PSO on a user specified datasegment and PSD. +%Refer to User Manual for usage + +%Raghav Girgaonkar, July 2023 + +%% Read JSON Files +fname = filename; +str = fileread(fname); +filenames = jsondecode(str); +fname = filenames.signalparamfile; +str = fileread(fname); +params = jsondecode(str); +fname = filenames.psoparamfile; +str = fileread(fname); +pso = jsondecode(str); +fname = filenames.filenames; +str = fileread(fname); +files = jsondecode(str); +%% Specify initial parameters +%T_sig_len = params.signal.T_sig_len; +T_sig = params.signal.T_sig; +initial_phase = 0; +%% Sampling Frequency +Fs = params.sampling_freq; +%% Number of samples = num*Fs*T_sig +%num = params.signal.num; +N = floor(Seglen*Fs); +timeVec = (0:N-1)*(1/Fs); +%% Min and Max Frequencies Hz +fmin = params.freq(1); +fmax = params.freq(2); +%% Positive Frequency Vector +datalen = N/Fs; +fpos = (0:floor(N/2))*(1/datalen); +%% Initial Time of Arrival and Coalescence Phase +ta = params.ta; +phase = params.phase; +%% Signal to noise ratio of the true signal +snr = params.snr; +%% Constants +c = 3*10^8; +Msolar = 1.989*10^30; +G = 6.6743*10^-11; +%% Mass parameters of the true signal +m1 = params.masses(1); +m2 = params.masses(2); +%% Tau coeffs as phase parameters +if pso.type == "tau" + m1 = m1*Msolar; + m2 = m2*Msolar; + M = m1 + m2; + u = m1*m2/(m1 + m2); + n = u/M; + tau0 = (5/(256*pi))*(1/fmin)*((G*M*pi*fmin/c^3)^(-5/3))*(1/n); + tau1p5 = (1/8)*(1/fmin)*((G*M*pi*fmin/c^3)^(-2/3))*(1/n); + type = 1; + disp("Tau Space PSO"); +else + type = 0; + disp("Mass Space PSO"); +end +%% Search range of phase coefficients +if type + rmin = [params.rmin_tau(1), params.rmin_tau(2)]; + rmax = [params.rmax_tau(1), params.rmax_tau(2)]; +else + rmin = [params.rmin(1), params.rmin(2)]; + rmax = [params.rmax(1), params.rmax(2)]; +end +%% Number of independent PSO runs +nRuns = pso.nruns; + +%% Pre-processing +[A,avec, phaseDiff] = preprocessing(fmin,fmax, fpos, datalen, N); + +dataY = whtnddata; + +AbysqrtPSD = A.*TFtotal; + +%% Create General Normalization Factor +% Scalar factor of 1/N (check NoteonNormalizaton.pdf for more details) +dataLen = N; +innProd = (1/dataLen)*(AbysqrtPSD)*AbysqrtPSD'; +genNormfacSqr = real(innProd); +genNormfac = 1/sqrt(genNormfacSqr); + +%% Data Products +fftdataY = fft(dataY); +fftdataY = fftdataY.*A;%Pre-multiply frequency magnitude vector A for optimization + +%% Get FFT of data by total PSD +fftdataYbyPSD = fftdataY.*TFtotal; + +dataX = timeVec; + +%% Input Parameters: +inParams = struct('dataX', dataX,... + 'fpos', fpos,... + 'dataY', dataY,... + 'fftdataYbyPSD', fftdataYbyPSD,... + 'frange', [fmin,fmax],... + 'datalen',datalen,..., + 'initial_phase', initial_phase,... + 'N', N,... + 'A', A,... + 'phaseDiff', phaseDiff,... + 'normfac', genNormfac,... + 'avec', avec,... + 'T_sig', T_sig,... + 'rmin',rmin,... + 'rmax',rmax,... + 'Fs',Fs); +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + +maxSteps = pso.maxSteps; +if type + if rmin(1) >= 0 && rmin(2) >= 0 + original_fitVal = -1*mfgw_tau([tau0, tau1p5], inParams); + outStruct = crcbgwpso_tau(inParams,struct('maxSteps',maxSteps),nRuns,Fs); + else + original_fitVal = -1*mfgw_tau_negative([tau0, tau1p5], inParams); + outStruct = crcbgwpso_tau_negative(inParams,struct('maxSteps',maxSteps),nRuns,Fs); + end + bestFitVal = -1*outStruct.bestFitness; +else + original_fitVal = -1*mfgw_mass([m1, m2], inParams); + outStruct = crcbgwpso_mass(inParams,struct('maxSteps',maxSteps),nRuns,Fs); + bestFitVal = -1*outStruct.bestFitness; +end + +%% Optional: Uncomment to save output of PSO in a .mat file +% save(files.output_struct_location,'outStruct'); + +%% Plots, uncomment all saveas() commands to save figures in custom directories +figure; +hold on; +plot(dataX,segdata,'.'); +% plot(dataX,wave,'r'); +% for lpruns = 1:nRuns +% plot(dataX,outStruct.allRunsOutput(lpruns).estSig,'Color',[51,255,153]/255,'LineWidth',4.0); +% end +plot(dataX,outStruct.bestSig,'Color',[76,153,0]/255,'LineWidth',2.0); +% legend('Data','Signal',... +% ['Estimated signal: ',num2str(nRuns),' runs'],... +% 'Estimated signal: Best run'); +legend('Data','Signal',... + 'Estimated signal: Best run'); +% saveas(gcf,files.psoresultplot); +hold off; + +figure; +iterVec = linspace(1,maxSteps,maxSteps); +hold on; +for lpruns = 1:nRuns + plot(iterVec,outStruct.allRunsOutput(lpruns).allBestFit, 'DisplayName',num2str(lpruns)); +end +title("Best Fitness Values for All Runs"); +xlabel("Iteration"); +ylabel("Best Fitness Value"); +legend; +% saveas(gcf,files.bestfitplot); +hold off; + +if type + figure; + hold on; + for lpruns = 1:nRuns + rVec = s2rv(outStruct.allRunsOutput(lpruns).allBestLoc,inParams); + plot(rVec(:,1),rVec(:,2),'DisplayName',num2str(lpruns)); + end + scatter(tau0,tau1p5,140,'red','filled','D','DisplayName','Original Parameters'); + title("Best Parameter Values for All Runs"); + xlabel("\tau_0"); + ylabel("\tau_{1.5}"); + legend; + boundary_plot; +% saveas(gcf,files.bestlocplot); + hold off; + + + + t0 = outStruct.bestQcCoefs(1); + t1p5 = outStruct.bestQcCoefs(2); + est_ta = outStruct.bestTime; + est_SNR = outStruct.bestAmp; + est_M = (5/(32*fmin))*(t1p5/(pi*pi*t0))*(c^3/G); + est_u = (1/(16*fmin*fmin))*(5/(4*pi^4*t0*t1p5^2))^(1/3)*(c^3/G); + + est_m1 = (est_M - sqrt(est_M^2 - 4*est_u*est_M))/(2*Msolar); + est_m2 = (est_M + sqrt(est_M^2 - 4*est_u*est_M))/(2*Msolar); + + %% This will display parameters given through signal.json and PSO-estimated parameters + %% Uncomment Original parameter display command if needed +% disp(['Original parameters: tau0= ',num2str(tau0),... +% '; tau1p5= ',num2str(tau1p5),... +% '; m1= ', num2str(m1/Msolar),... +% '; m2= ', num2str(m2/Msolar),... +% '; A = ',num2str(snr),... +% '; phi = ',num2str(phase),... +% '; t_a = ',num2str(ta),... +% '; FitVal = ',num2str(original_fitVal)]); + + disp(['Estimated parameters: tau0=',num2str(outStruct.bestQcCoefs(1)),... + '; tau1p5=',num2str(outStruct.bestQcCoefs(2)),... + '; m1= ', num2str(est_m1),... + '; m2= ', num2str(est_m2),... + '; A = ',num2str(outStruct.bestAmp),... + '; phi = ',num2str(outStruct.bestPhase),... + '; t_a = ',num2str(outStruct.bestTime),... + '; FitVal = ',num2str(bestFitVal)]); +else + figure; + hold on; + for lpruns = 1:nRuns + rVec = s2rv(outStruct.allRunsOutput(lpruns).allBestLoc,inParams); + plot(rVec(:,1),rVec(:,2),'DisplayName',num2str(lpruns)); + end + scatter(m1,m2,140,'red','filled','D','DisplayName','Original Parameters'); + title("Best Parameter Values for All Runs"); + xlabel("m_1"); + ylabel("m_2"); + legend; +% saveas(gcf,files.bestlocplot); + hold off; + + %% This will display parameters given through signal.json and PSO-estimated parameters + %% Uncomment Original parameter display command if needed +% disp(['Original parameters: m1= ',num2str(m1),... +% '; m2= ',num2str(m2),... +% '; A = ',num2str(snr),... +% '; phi = ',num2str(phase),... +% '; t_a = ',num2str(ta),... +% '; FitVal = ',num2str(original_fitVal)]); + + est_m1 = outStruct.bestQcCoefs(1); + est_m2 = outStruct.bestQcCoefs(2); + est_ta = outStruct.bestTime; + est_SNR = outStruct.bestAmp; + disp(['Estimated parameters: m1=',num2str(outStruct.bestQcCoefs(1)),... + '; m2=',num2str(outStruct.bestQcCoefs(2)),... + '; A = ',num2str(outStruct.bestAmp),... + '; phi = ',num2str(outStruct.bestPhase),... + '; t_a = ',num2str(outStruct.bestTime),... + '; FitVal = ',num2str(bestFitVal)]); +end + +%{ +S = load('injectedsignals.mat'); sigsegnum = S.segnum; sigsnrs = S.injectedsnrs; + +index = find(sigsegnum == segnum); + +if index + signalsnr = sigsnrs(index); +else + signalsnr = 0; +end +%} + + +%candidates = [segnum, t0, t1p5, real(est_m1), imag(est_m1), real(est_m2), imag(est_m2), est_ta, est_SNR,signalsnr]; +candidates = [segnum, t0, t1p5, real(est_m1), imag(est_m1), real(est_m2), imag(est_m2), est_ta, est_SNR]; +%candidates = [segnum, est_m1, est_m2, est_ta, est_SNR]; +fileID = fopen(files.candidatefile,'a'); +fprintf(fileID,'%d\t%f\t%f\t%f + %fi\t\t%f + %fi\t\t%f\t\t%f\n',candidates); +%fprintf(fileID,'%d\t%f\t%f\t%f + %fi\t\t%f + %fi\t\t%f\t\t%f\t%f\n',candidates); +%fprintf(fileID,'%d\t%f\t%f\t%f\t%f\n',candidates); +fclose(fileID); +end + + diff --git a/OptimizedPSOCodes/segdatacond.m b/OptimizedPSOCodes/segdatacond.m index 02ea78b..f8a66e6 100644 --- a/OptimizedPSOCodes/segdatacond.m +++ b/OptimizedPSOCodes/segdatacond.m @@ -1,4 +1,4 @@ -function [whtndfiltdata, whtndstd, TFtotal]=segdatacond(segdata, PSD, sampFreq, tidxs) +function [whtndfiltdata, TFtotal]=segdatacond(segdata, PSD, sampFreq) fmin = 30; % sampfreq = 4096; @@ -38,9 +38,12 @@ whtndfiltdata = ifft(whtndfftfiltdata); -%Divide by variance of whitened strain so that final whitened vector has -%unit strain -tstart = tidxs(1); -tend = tidxs(2); -whtndstd = std(whtndfiltdata(tstart: tend)); -whtndfiltdata = whtndfiltdata/whtndstd; +%Normalize it to unit variance +whtndfiltdata = whtndfiltdata/sqrt(sampFreq); + +% %Divide by variance of whitened strain so that final whitened vector has +% %unit strain +% tstart = tidxs(1); +% tend = tidxs(2); +% whtndstd = std(whtndfiltdata(tstart: tend)); +% whtndfiltdata = whtndfiltdata/whtndstd; diff --git a/animate_matchedfiltering.m b/animate_matchedfiltering.m new file mode 100644 index 0000000..4d8ce45 --- /dev/null +++ b/animate_matchedfiltering.m @@ -0,0 +1,127 @@ +%% Script to create an animation of matched filtering timeseries. +addpath('Misc/NormalizationCheck/') +%% Define Data Parameters +sampFreq = 4096; +datalen = 128; +N = datalen*sampFreq; +frange = [30,700]; +fpos = (0:floor(N/2))*(1/datalen); +timeVec = (0:(datalen*sampFreq - 1))/sampFreq; + +%% Create Colored Gaussian Noise +[noise, PSD] = LIGOnoise(N, sampFreq, 1, 'sample'); + + +%Create Total PSD and Transfer Function +negFStrt = 1-mod(N,2); +kNyq = floor(N/2)+1; + +PSDtotal = [PSD, PSD((kNyq-negFStrt):-1:2)]; +TFtotal = 1./sqrt(PSDtotal); + + +%Pre-calculate frequency magnitude and phase difference terms + [A,avec, phaseDiff] = preprocessing(frange(1),frange(2),fpos, datalen, N); + +%Create general normalization factor N_fW +AbysqrtPSD = A.*TFtotal; +innProd = (1/N)*(AbysqrtPSD)*AbysqrtPSD'; +genNormfacSqr = real(innProd); +N_fW = 1/sqrt(genNormfacSqr); + +%% Generate Signal to Inject + chirptimes = [29.6373, 1.1045]; + ta = 34; + initial_phase = 0; + phase = 0; + snr = 12; + +signal = gen2PNtemplate_tau(fpos, ta, phase, frange(1), frange(2),chirptimes(1), chirptimes(2),datalen,initial_phase,snr, N, A, avec, PSDtotal); + +signal = signal*sqrt(sampFreq); + +figure; +plot(timeVec, noise); +hold on; +plot(timeVec, 30*signal); +xlabel('Time (s)'); +ylabel('Strain h(t)'); +axis tight; +ax = gca(); ax.XAxis.FontSize = 20; ax.YAxis.FontSize = 20; +%Save whitened and normalized signal and noise +%Window first +whtndsignal = ifft((1/sqrt(sampFreq))*fft(signal).*TFtotal); + +%Inject signal into data +strain = noise + signal; + +%Whiten and normalize strain +rolloff = 0.5; +fmin = frange(1); +strainwin = strain.*tukeywin(length(strain),rolloff*sampFreq/length(strain))'; +whtndstrain = ifft((1/sqrt(sampFreq))*fft(strainwin).*TFtotal); + + +%Make GIF +% Define the filename for the gif +giffilename = 'matchedfiltering_animation.gif'; +time_lapse = 0.0001; +title_text = 'Generation of Matched-Filtering Timeseries'; + +% Create a new figure window +fig = figure; +mftimeseries = zeros(1,N); +for i = 1:1*sampFreq:N + q0 = gen2PNtemplate_tau(fpos, i/sampFreq, 0, frange(1), frange(2),chirptimes(1), chirptimes(2),datalen,initial_phase,1, N, A, avec, PSDtotal); + whtndq0 = ifft(fft(q0).*TFtotal); + q1 = gen2PNtemplate_tau(fpos, i/sampFreq, pi/2, frange(1), frange(2),chirptimes(1), chirptimes(2),datalen,initial_phase,1, N, A, avec, PSDtotal); + whtndq1 = ifft(fft(q1).*TFtotal); + +% I0 = innerproduct(strain, q0, PSDtotal, sampFreq); +% I1 = innerproduct(strain, q1, PSDtotal, sampFreq); + + I0 = real((1/N)*sum(fft(whtndstrain).*conj(fft(whtndq0)))); + I1 = real((1/N)*sum(fft(whtndstrain).*conj(fft(whtndq1)))); + + I = sqrt(I0^2 + I1^2); + + mftimeseries(i) = I; + + + hold on; + plot(timeVec, whtndstrain); hold on; + plot(timeVec, 5*whtndsignal, Color='red' ); hold on; + plot(timeVec, 20*whtndq0,Color='yellow'); hold on; + plot(timeVec(1:i), mftimeseries(1:i), Color='black',LineWidth=1.5); hold on; + ylim([-4,12]); + xlim([0,128]); + xlabel("Time (s)"); + ylabel("Whitened Strain h(t)"); + ax = gca(); ax.XAxis.FontSize = 20; ax.YAxis.FontSize = 20; +% axis tight; + + title(title_text,'FontSize',16,'Color',[0 0 0]), hold on; + + drawnow; + %Capture frame + % Capture the current frame and save it to the gif + frame = getframe(fig); + im = frame2im(frame); + [imind,cm] = rgb2ind(im,256); + + if i == 1 + imwrite(imind,cm,giffilename,'gif', 'Loopcount',Inf,'DelayTime',time_lapse); + else + imwrite(imind,cm,giffilename,'gif','WriteMode','append','DelayTime',time_lapse); + end + + clf(fig); + +end + +% +% figure; +% plot(timeVec, noise); +% hold on; +% plot(timeVec, 10*signal); +% axis tight; diff --git a/animate_matchedfiltering_glitch.m b/animate_matchedfiltering_glitch.m new file mode 100644 index 0000000..0f5fa9f --- /dev/null +++ b/animate_matchedfiltering_glitch.m @@ -0,0 +1,128 @@ +%% Script to create an animation of matched filtering timeseries. +addpath('Misc/NormalizationCheck/') +%% Define Data Parameters +sampFreq = 4096; +datalen = 128; +N = datalen*sampFreq; +frange = [30,700]; +fpos = (0:floor(N/2))*(1/datalen); +timeVec = (0:(datalen*sampFreq - 1))/sampFreq; + +%% Create Colored Gaussian Noise +[noise, PSD] = LIGOnoise(N, sampFreq, 1, 'sample'); + + +%Create Total PSD and Transfer Function +negFStrt = 1-mod(N,2); +kNyq = floor(N/2)+1; + +PSDtotal = [PSD, PSD((kNyq-negFStrt):-1:2)]; +TFtotal = 1./sqrt(PSDtotal); + + +%Pre-calculate frequency magnitude and phase difference terms + [A,avec, phaseDiff] = preprocessing(frange(1),frange(2),fpos, datalen, N); + +%Create general normalization factor N_fW +AbysqrtPSD = A.*TFtotal; +innProd = (1/N)*(AbysqrtPSD)*AbysqrtPSD'; +genNormfacSqr = real(innProd); +N_fW = 1/sqrt(genNormfacSqr); + +%% Generate Signal to Inject + chirptimes = [0.001,0.01]; + qchirptimes = [11.6373, 1.1045]; + ta = 54; + initial_phase = 0; + phase = 0; + snr = 360; + +signal = gen2PNtemplate_tau_negative(fpos, ta, phase, frange(1), frange(2),chirptimes(1), chirptimes(2),datalen,initial_phase,snr, N, A, avec, PSDtotal); + +signal = signal*sqrt(sampFreq); + +% figure; +% plot(timeVec, noise); +% hold on; +% plot(timeVec, 30*signal); +% xlabel('Time (s)'); +% ylabel('Strain h(t)'); +% axis tight; +% ax = gca(); ax.XAxis.FontSize = 20; ax.YAxis.FontSize = 20; +%Save whitened and normalized signal and noise +%Window first +whtndsignal = ifft((1/sqrt(sampFreq))*fft(signal).*TFtotal); + +%Inject signal into data +strain = noise + signal; + +%Whiten and normalize strain +rolloff = 0.5; +fmin = frange(1); +strainwin = strain.*tukeywin(length(strain),rolloff*sampFreq/length(strain))'; +whtndstrain = ifft((1/sqrt(sampFreq))*fft(strainwin).*TFtotal); + + +%Make GIF +% Define the filename for the gif +giffilename = 'matchedfiltering_glitch_animation.gif'; +time_lapse = 0.0001; +title_text = 'Generation of Matched-Filtering Timeseries'; + +% Create a new figure window +fig = figure; +mftimeseries = zeros(1,N); +for i = 1:1*sampFreq:N + q0 = gen2PNtemplate_tau(fpos, i/sampFreq, 0, frange(1), frange(2),qchirptimes(1), qchirptimes(2),datalen,initial_phase,1, N, A, avec, PSDtotal); + whtndq0 = ifft(fft(q0).*TFtotal); + q1 = gen2PNtemplate_tau(fpos, i/sampFreq, pi/2, frange(1), frange(2),qchirptimes(1), qchirptimes(2),datalen,initial_phase,1, N, A, avec, PSDtotal); + whtndq1 = ifft(fft(q1).*TFtotal); + +% I0 = innerproduct(strain, q0, PSDtotal, sampFreq); +% I1 = innerproduct(strain, q1, PSDtotal, sampFreq); + + I0 = real((1/N)*sum(fft(whtndstrain).*conj(fft(whtndq0)))); + I1 = real((1/N)*sum(fft(whtndstrain).*conj(fft(whtndq1)))); + + I = sqrt(I0^2 + I1^2); + + mftimeseries(i) = I; + + + hold on; + plot(timeVec, whtndstrain); hold on; + plot(timeVec, whtndsignal, Color='red' ); hold on; + plot(timeVec, 20*whtndq0,Color='yellow'); hold on; + plot(timeVec(1:i), mftimeseries(1:i), Color='black',LineWidth=1.5); hold on; + ylim([-4,15]); + xlim([0,128]); + xlabel("Time (s)"); + ylabel("Whitened Strain h(t)"); + ax = gca(); ax.XAxis.FontSize = 20; ax.YAxis.FontSize = 20; +% axis tight; + + title(title_text,'FontSize',16,'Color',[0 0 0]), hold on; + + drawnow; + %Capture frame + % Capture the current frame and save it to the gif + frame = getframe(fig); + im = frame2im(frame); + [imind,cm] = rgb2ind(im,256); + + if i == 1 + imwrite(imind,cm,giffilename,'gif', 'Loopcount',Inf,'DelayTime',time_lapse); + else + imwrite(imind,cm,giffilename,'gif','WriteMode','append','DelayTime',time_lapse); + end + + clf(fig); + +end + +% +% figure; +% plot(timeVec, noise); +% hold on; +% plot(timeVec, 10*signal); +% axis tight;