Added negative chirp waveforms

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)");

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);
 

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}');

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;
+

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;
+PSD = (10.^loginterPSD);
+% PSD = (loginterPSD)./2;

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;
+
+
+

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;

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
+

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;