P = cymbal_defaults(); Q = P; Q.normalize = 0; Q.full_scale_pa = 1;
fs = P.sample_rate; win = 4096;
for turb = [0 0.9]
  R = Q; R.turbulence = turb;
  [y, b, d] = cymbal_render(R);
  fprintf('turbulence %.1f  peak %6.2f Pa  flat@ ', turb, d.peakPa);
  for tt = [0.01 0.10 0.40 1.00 2.00]
    i0 = round(tt*fs)+1;
    x = y(i0:i0+win-1) .* (0.5*(1-cos(2*pi*(0:win-1)/(win-1))));
    Y = abs(fft(x)).^2; Y = Y(1:win/2); fq = (0:win/2-1)*fs/win;
    Ys = Y(fq>=1000 & fq<8000) + 1e-24;
    fprintf('%.3f ', exp(mean(log(Ys)))/mean(Ys));
  end
  i0 = round(1.0*fs)+1;
  x = y(i0:i0+win-1) .* (0.5*(1-cos(2*pi*(0:win-1)/(win-1))));
  Y = abs(fft(x)).^2; Y = Y(1:win/2); fq = (0:win/2-1)*fs/win;
  fprintf(' | 1s low-vs-high %.0f dB\n', 10*log10(sum(Y(fq<600))+1e-30) - 10*log10(sum(Y(fq>=1500&fq<8000))+1e-30));
end
P2 = cymbal_defaults();
[y2, b2, d2] = cymbal_render(P2);
fprintf('\npinned: peakPa %.6f  f1 %.6f  n %d\n', d2.peakPa, b2.f(1), b2.n);