測量鋸齒的週期

Question

我有一系列二維測量結果（x軸上的時間），繪製成非光滑（但非常好）的鋸齒波。在理想的世界中，數據點會形成完美的鋸齒波（在任一端都有部分振幅數據點）。有沒有使用OCTAVE/MATLAB來計算波的（平均）週期的方法？我嘗試使用公式維基百科（Sawtooth_wave）鋸齒：

P = mean(time.*pi./acot(tan(y./4))), -pi < y < +pi 

也試過：

P = mean(abs(time.*pi./acot(tan(y./4)))) 

，但它沒有工作，或至少它給了我，我知道的答案了。

繪製數據的一個例子：

我也嘗試下面的方法 - 應該工作 - 但它不給我我所知道的是接近正確的答案。可能是我的代碼簡單和錯誤。什麼？

slopes = diff(y)./diff(x); % form vector of slopes for each two adjacent points 
for n = 1:length(diff(y)) % delete slope of any two points that form the 'cliff' 
    if abs(diff(y(n,1))) > pi 
    slopes(n,:) = []; 
    end 
    end 
P = median((2*pi)./slopes); % Amplitude is 2*pi 

Answer 1

舊帖子，但認爲我會提供我的2美分的價值。我認爲有兩個合理的方式來做到這一點：

1. 進行傅立葉變換，計算基本
2. 執行階段，週期，幅度的曲線擬合，並且偏移到理想的方波。

鑑於曲線擬合可能很困難，因爲鋸齒波的不連續性，所以我建議傅里葉變換。下面是一個獨立的示例：

f_s = 10;    # Sampling freq. in Hz 
record_length = 1000; # length of recording in sec. 

% Create noisy saw-tooth wave, with known period and phase 
saw_period = 50; 
saw_phase = 10; 
t = (1/f_s):(1/f_s):record_length; 
saw_function = @(t) mod((t-saw_phase)*(2*pi/saw_period), 2*pi) - pi; 

noise_lvl = 2.0; 
saw_wave = saw_function(t) + noise_lvl*randn(size(t)); 
num_tsteps = length(t); 

% Plot time-series data 
figure(); 
plot(t, saw_wave, '*r', t, saw_function(t)); 
xlabel('Time [s]'); 
ylabel('Measurement'); 
legend('measurements', 'ideal'); 

% Perform fast-Fourier transform (and plot it) 
dft = fft(saw_wave); 
freq = 0:(f_s/length(saw_wave)):(f_s/2); 
dft = dft(1:(length(saw_wave)/2+1)); 

figure(); 
plot(freq, abs(dft)); 
xlabel('Freqency [Hz]'); 
ylabel('FFT of Measurement'); 

% Estimate fundamental frequency: 
[~, idx] = max(abs(dft)); 
peak_f = abs(freq(idx)); 
peak_period = 1/peak_f; 
disp(strcat('Estimated period [s]: ', num2str(peak_period))) 

其中輸出一對圖表，以及鋸齒波的估計週期。你可以玩弄噪音的數量，並看到它正確地獲得了50秒的時間，直到非常高的噪音水平。

Estimated period [s]: 50