From b4a8b2bd106e95d2875517efd77558867368a4a0 Mon Sep 17 00:00:00 2001 From: RaghavGirgaonkar Date: Thu, 14 Dec 2023 15:52:05 -0600 Subject: [PATCH] Added negative chirp waveforms --- .../AllPSO/animate_matchedfiltering.m | 4 +- OptimizedPSOCodes/LIGOnoise.m | 4 +- OptimizedPSOCodes/boundary_plot.m | 12 +- OptimizedPSOCodes/crcbgwpso_tau_negative.m | 123 ++++++++++++++++++ OptimizedPSOCodes/createPSD.m | 42 ++---- .../gen2PNwaveform_tau_negative.m | 20 +++ OptimizedPSOCodes/mfgw_tau_negative.m | 27 ++++ OptimizedPSOCodes/psofitfunc_tau_negative.m | 38 ++++++ OptimizedPSOCodes/waveform_tau.m | 6 +- OptimizedPSOCodes/waveform_tau_negative.m | 102 +++++++++++++++ 10 files changed, 340 insertions(+), 38 deletions(-) create mode 100644 OptimizedPSOCodes/crcbgwpso_tau_negative.m create mode 100644 OptimizedPSOCodes/gen2PNwaveform_tau_negative.m create mode 100644 OptimizedPSOCodes/mfgw_tau_negative.m create mode 100644 OptimizedPSOCodes/psofitfunc_tau_negative.m create mode 100644 OptimizedPSOCodes/waveform_tau_negative.m diff --git a/Misc/2PNWaveform/AllPSO/animate_matchedfiltering.m b/Misc/2PNWaveform/AllPSO/animate_matchedfiltering.m index 11d0226..6d935a7 100644 --- a/Misc/2PNWaveform/AllPSO/animate_matchedfiltering.m +++ b/Misc/2PNWaveform/AllPSO/animate_matchedfiltering.m @@ -1,6 +1,6 @@ function []=animate_matchedfiltering(filename) %% Read JSON File -% addpath("../../PSO/"); +addpath("../../AllPSO/"); fname = filename; str = fileread(fname); filenames = jsondecode(str); @@ -137,7 +137,7 @@ hold on; plot(timeVec, dataY); hold on; plot(timeVec, wave, Color='red' ); hold on; -plot(timeVec, q0,Color='yellow'); hold on; +plot(timeVec, 3*q0,Color='yellow'); hold on; plot(timeVec(1:i), mftimeseries(1:i), Color='black'); hold on; ylim([-4,12]); xlabel("Time (s)"); diff --git a/OptimizedPSOCodes/LIGOnoise.m b/OptimizedPSOCodes/LIGOnoise.m index d4c242c..66d9b45 100644 --- a/OptimizedPSOCodes/LIGOnoise.m +++ b/OptimizedPSOCodes/LIGOnoise.m @@ -39,7 +39,7 @@ interPSD = interp1(y(:,1),y(:,2), fvec); %% Modifications, change cutoff frequencies as needed -minidx = find(fvec<=15, 1, 'last' ); +minidx = find(fvec<=30, 1, 'last' ); maxidx = find(fvec<=700, 1, 'last' ); Sn50 = interPSD(minidx); @@ -53,7 +53,7 @@ %% Make colored Noise fltrOrdr = 10000; -outNoise_t = statgaussnoisegen(N+ 2*fltrOrdr,[fvec(:),PSD(:)],fltrOrdr,Fs, noise_num, noisefile); +outNoise_t = statgaussnoisegen(N,[fvec(:),PSD(:)],fltrOrdr,Fs, noise_num, noisefile); outNoise = outNoise_t(fltrOrdr+1:end - fltrOrdr); diff --git a/OptimizedPSOCodes/boundary_plot.m b/OptimizedPSOCodes/boundary_plot.m index 78ce0e2..b268220 100644 --- a/OptimizedPSOCodes/boundary_plot.m +++ b/OptimizedPSOCodes/boundary_plot.m @@ -1,4 +1,4 @@ -function boundary_fig = boundary_plot +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 @@ -19,7 +19,7 @@ M = (m1 + m2); u = m1.*m2./M; n = u./M; -%Calculate Chirp Times +%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)); @@ -63,6 +63,14 @@ %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}'); diff --git a/OptimizedPSOCodes/crcbgwpso_tau_negative.m b/OptimizedPSOCodes/crcbgwpso_tau_negative.m new file mode 100644 index 0000000..2e4c68a --- /dev/null +++ b/OptimizedPSOCodes/crcbgwpso_tau_negative.m @@ -0,0 +1,123 @@ +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); + [~,qcCoefs,ta_index] = fHandle(outStruct(lpruns).bestLocation); +% Q = qcCoefs +% index = ta_index + outResults.allRunsOutput(lpruns).qcCoefs = qcCoefs; + %Calculate time using sampling freq and ta_index + estTa = ta_index/sampling_freq; + + outResults.allRunsOutput(lpruns).estTa = estTa; + tau0 = qcCoefs(1); + tau1p5 = qcCoefs(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/OptimizedPSOCodes/createPSD.m b/OptimizedPSOCodes/createPSD.m index 5b0583a..60bd652 100644 --- a/OptimizedPSOCodes/createPSD.m +++ b/OptimizedPSOCodes/createPSD.m @@ -1,43 +1,23 @@ -%Script to create interpolated PSD from SHAPES and PWELCH Estimates - -%% Provide Filename containing log-estimates for SHAPES and PWELCH -S = load('~/Documents/linetest_Job_col_PSD_dataSz4Up2N1024SNRR5e-03_C.mat'); - -%% Get One-sided estimates - -welchPSD = S.dataY; -shapesPSD = S.estS; +function [PSD, fvec]=createPSD(sampFreq, Tsig, welchPSD, freqs) +%% Script to create interpolated PSD from SHAPES and PWELCH Estimates %% Data Parameters -Fs = 4096; -T = 4096; +Fs = sampFreq; +T = Tsig; N = Fs*T; kNyq = floor(N/2); fvec = (0:(kNyq))*Fs/N; -minidx = find(fvec<=15, 1, 'last' ); -maxidx = find(fvec<=700, 1, 'last' ); - -nsamp = kNyq - minidx + 1; +% nyqfreq = winlen*Fs; -%% Complete the PSD Vectors - -winlen = 4; %Window length of Welch Estimate -n = winlen*Fs; - -shapespsd_temp = [ones(1, n+1-length(shapesPSD))*shapesPSD(1), shapesPSD']; -welchpsd_temp = [ones(1, n+1-length(shapesPSD))*welchPSD(1), welchPSD']; - -nyqfreq = n; -freqs = (0:nyqfreq)*(Fs/(2*n)); +% freqs = (0:nyqfreq)*(Fs/(2*nyqfreq)); %% 1-D Interpolation +logwelchPSD = log10(welchPSD); +loginterPSD = interp1(freqs, logwelchPSD, fvec); -WPSD = interp1(freqs, welchpsd_temp, fvec); -SPSD = interp1(freqs, shapespsd_temp, fvec); - -%% Antilog +% %% Antilog -SHAPESPSD = (10.^SPSD)./2; -WELCHPSD = (10.^WPSD)./2; \ No newline at end of file +PSD = (10.^loginterPSD); +% PSD = (loginterPSD)./2; diff --git a/OptimizedPSOCodes/gen2PNwaveform_tau_negative.m b/OptimizedPSOCodes/gen2PNwaveform_tau_negative.m new file mode 100644 index 0000000..f3f54c4 --- /dev/null +++ b/OptimizedPSOCodes/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/OptimizedPSOCodes/mfgw_tau_negative.m b/OptimizedPSOCodes/mfgw_tau_negative.m new file mode 100644 index 0000000..ffeaf06 --- /dev/null +++ b/OptimizedPSOCodes/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/OptimizedPSOCodes/psofitfunc_tau_negative.m b/OptimizedPSOCodes/psofitfunc_tau_negative.m new file mode 100644 index 0000000..3074f6e --- /dev/null +++ b/OptimizedPSOCodes/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/OptimizedPSOCodes/waveform_tau.m b/OptimizedPSOCodes/waveform_tau.m index 10601bd..2dc9b34 100644 --- a/OptimizedPSOCodes/waveform_tau.m +++ b/OptimizedPSOCodes/waveform_tau.m @@ -17,7 +17,11 @@ G = 6.6743*10^-11; %% Calculate Mass Terms from Tau0 and Tau1.5 M = (5/(32*fmin))*(tau1p5/(pi*pi*tau0))*(c^3/G); -u = (1/(16*fmin*fmin))*(5/(4*pi^4*tau0*tau1p5^2))^(1/3)*(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; diff --git a/OptimizedPSOCodes/waveform_tau_negative.m b/OptimizedPSOCodes/waveform_tau_negative.m new file mode 100644 index 0000000..db53603 --- /dev/null +++ b/OptimizedPSOCodes/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 +if abs(tau0) < abs(tau1p5) + 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 +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; + + + + +