Here are two simple approaches that utilize all the samples to determine the result and work under much noisier conditions with a minimum error in the result. Both of these methods are significantly more robust under conditions of noise than any edge detection methods as they are using all samples, not just the samples in proximity to an edge.
First we need to minimize the error due to truncating the square wave at a non-integer boundary of the cycle rate. If the square wave is of high fidelity as shown in the OP's graph, this can be done by simply truncating the ending samples such that the transition from the very last sample in the graph is continuous to the first sample (if we put the graph end on end, there would be no discontinuity). If that isn't convenient, the other approach is to window the waveform so that the end points where the discontinuity exists do not contribute as much to the result, and we concentrate our efforts in the middle of the graph.
For example, using wfm
to represent the square wave:
wfm = wfm(:); # force to be column vector
N = length(wfm);
win = kaiser(N, 8);
wfm_win = wfm .* win;
wfm_win
is used in each approach below to measure the phase between the square wave and a sinusoid starting at phase = 0 at t=0.
FFT Approach
Zero pad wfm_win
to interpolate the time domain result, and then select the maximum bin of the resulting FFT. The phase of this bin is the desired phase measurement. The zero padding will minimize offset errors between the frequency of the maximum bin and the actual exact frequency of the square wave.
Phase Detector Approach
If you already know the exact frequency of the waveform and just want to get the phase, multiply the waveform by samples of $e^{-j 2 \pi f t}$ where f is the frequency. Take the average of the result (to low pass filter out the higher frequency) and then use atan2(Q,I)
to get the phase of the result, where $I$ and $Q$ are the real and imaginary terms of the multiplier output respectively. The first approach and the second approach are identical, to the extent that the interpolated bin in the FFT matches the "exact frequency" as used in the second approach since this here would be the computation of one FFT bin.
Demonstration
Here is a Matlab demonstration of both approaches uses a waveform that visually appears similar to the OP's waveform, with a square wave starting at $t=0$ and thus the resulting expected phase for a sine wave would be 90°.
t = linespace(0,47, 1000);
wfm = square(t)';
win = kaiser(length(wfm), 8);
wfm_win = wfm.*win;
# method 1
out = fft(wfm_win, 2**18); # zero pads wfm_win out to 2**18 samples
[val, idx]= max(out);
ph_result = atan2(imag(out(idx)), real(out(idx)))
# method 2
ref = exp(-1j*t)';
out2 = wfm_win.*ref;
ph_result2 = atan2(imag(mean(out2)), real(mean(out2)))
The result in this case was 1.5698 radians (89.944 degrees) for method 1 and 1.5683 radians (89.857 degrees) for method 2.
Plot of windowed square wave used:

T
, which may very well produce some incorrect (or not expected) result which propagates topulsewidth()
function and produces an error. Don't take my word for it though, this is just an idea, since I have never usedpulsewidth()
before. $\endgroup$.wav
. There are also tricks in forcing the timestep to be sharp around the samples of interest, using aPULSE()
source (similar to what you have there). But, if you want a sampled waveform then why not use Matlab/Octave, directly? $\endgroup$