From 921fec4c54c0ca91f3caf86aae60385f3e89a20a Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Mon, 31 Jan 2022 22:20:24 +0100 Subject: [PATCH 1/8] lag_energy_spectrum: Improve run-time for freq-resolved LES * Previously, frequency-resolved lag-energy spectrum could be completed by looping over the frequency bins and then calling lag_energy_spectrum, which loops through the energy bins. However, it is much faster to first loop through the energies and then through the frequencies (the DFT has to be done only once for each LC). * Change f_min, f_max qualifiers to be arrays (containing the frequency grid) instead of doubles. The loop through the frequency bins is now implemented in the lag_energy_spectrum function * Change output struct: now contains all available information. This includes a change of the field errlag -> errlag_nowak99 * Export file writing into own function, remove option to pass open file pointer to lag_energy_spectrum function --- src/data/timing/lag_energy_spectrum.sl | 267 +++++++++++++------------ 1 file changed, 140 insertions(+), 127 deletions(-) diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl index 79fef22a..cdf35b40 100644 --- a/src/data/timing/lag_energy_spectrum.sl +++ b/src/data/timing/lag_energy_spectrum.sl @@ -3,11 +3,12 @@ private define filterFouQuantForFrequencyRange(fouquan, f_min, f_max) { %{{{ variable filt = where(f_min < fouquan.freq <= f_max); - if (qualifier("verbose", 0) > 1){ - vmessage(" ***(lag_energy_spectrum): Minimal frequency bin %d at %g Hz", filt[0], fouquan.freq[filt[0]]); - vmessage(" ***(lag_energy_spectrum): Maximal frequency bin %d at %g Hz", filt[-1], fouquan.freq[filt[-1]]); - vmessage(" ***(lag_energy_spectrum): Number of frequency bins: %d", length(filt)); - } + if (length(filt) == 0) vmessage("ERROR (lag_energy_spectrum): No frequency bins between %g-%g in Fourier products!", f_min, f_max); + if (qualifier("verbose", 0) > 1){ + vmessage(" ***(lag_energy_spectrum): Minimal frequency bin %d at %g Hz", filt[0], fouquan.freq[filt[0]]); + vmessage(" ***(lag_energy_spectrum): Maximal frequency bin %d at %g Hz", filt[-1], fouquan.freq[filt[-1]]); + vmessage(" ***(lag_energy_spectrum): Number of frequency bins: %d", length(filt)); + } return struct_filter(fouquan, filt ; copy); } %}}} @@ -89,6 +90,61 @@ private define write_cpd_to_ascii(realcpd, realcpd_err, imagcpd, imagcpd_err, f_ () = fclose(fptr_real); } %}}} +private define write_to_outfile(dat, f_lo, f_hi, dimseg, dt, normtype, lc_list){ + %{{{ + % To Do: MIN_FREQS, MAX_FREQS, N_FREQUENCIES, and N_SEGMENTS should + % be written as keywords, not columns, because they're the same + % within each frequency extension. + + variable fp = fits_open_file(qualifier("outfile"), "c"); + + variable N_energies, N_frequencies; + (N_energies, N_frequencies) = array_shape(dat.lag)[0], array_shape(dat.lag)[1]; + + variable jj; + _for jj (0, N_frequencies-1, 1){ + % write_cpd_to_ascii(dat.mean_realcpds[*,jj], dat.mean_realcpd_errs[*,jj], + % dat.mean_imagcpds[*,jj], dat.mean_imagcpd_errs[*,jj], + % f_lo[jj], f_hi[jj]); + + variable extname = sprintf("BAND_%f_%fHz", f_lo[jj], f_hi[jj]); + variable sdata = struct{ + "ENERGY_LO" = qualifier("elo", Double_Type[N_energies]), + "ENERGY_HI" = qualifier("ehi", Double_Type[N_energies]), + "LAG" = dat.lag[*,jj], + "LAG_ERR_UTTLEY14" = dat.errlag_uttley14[*,jj], + "LAG_ERR_NOWAK99" = dat.errlag_nowak99[*,jj], + "LAG_ERR_PROPERR" = dat.errlag_properr[*,jj], + "LAG_ERR_INGRAM19" = dat.errlag_ingram19[*,jj], + "MEAN_IMAG_CPD" = dat.mean_imagcpds[*,jj], + "MEAN_IMAG_CPD_ERR" = dat.mean_imagcpd_errs[*,jj], + "MEAN_REAL_CPD" = dat.mean_realcpds[*,jj], + "MEAN_REAL_CPD_ERR" = dat.mean_realcpd_errs[*,jj], + "MEAN_PHASE" = dat.mean_phases[*,jj], + "N_SEGMENTS" = dat.n_segments[*,jj], + "N_FREQUENCIES" = dat.n_frequencies[*,jj], + "MIN_FREQUENCY" = dat.min_freqs[*,jj], + "MAX_FREQUENCY" = dat.max_freqs[*,jj], + "RMS" = dat.rms, + "RMS_REF" = dat.rms_refs + }; + variable keys = struct{ + F_MIN=f_lo[jj], F_MAX=f_hi[jj], DIMSEG=dimseg, TIMEDEL=dt, NORMTYPE=normtype + }; + variable hist = struct{ + history = [sprintf("Wrote lag-energy spectrum in frequency range %f-%fHz", f_lo[jj], f_hi[jj])], + comment = ["Time unit is seconds, energy unit is keV, frequency unit is Hertz", + "Analyzed lightcurves:", lc_list] + }; + fits_write_binary_table(fp, extname, sdata, keys, hist); + } + + fits_close_file(fp); + () = system(sprintf("fthedit %s\[0] keyword=N_FREQS operation=add value=%i comment=\"Number of frequency bands\"", + qualifier("outfile"), N_frequencies)); + vmessage("Wrote lag-energy spectrum data to %s", qualifier("outfile")); +} %}}} + define lag_energy_spectrum(lc_list, dimseg) { %!%+ @@ -110,14 +166,13 @@ define lag_energy_spectrum(lc_list, dimseg) { % dt - Time resolution of the lightcurve in seconds % [default: TIMEDEL keyword] % deadtime - Detector deadtime in seconds [default: 0.0s] -% f_min - For frequency-resolved. Default: lowest sampled frequency -% [1/(dt*dimseg)]. Unit: Hz -% f_max - See f_min. Default: Nyquist frequency [1/(2*dt)]. Unit: Hz +% f_lo - Array_Type. For frequency-resolved lag-energy +% spectra. Default: lowest sampled frequency [an array +% containing only f_min=1/(dt*dimseg)]. Unit: Hz +% f_hi - See f_lo. Default: Nyquist frequency [1/(2*dt)]. Unit: Hz % verbose - Increase verbosity [default=0] -% outfile - String_Type or Fits_File_Type. If set, write output -% either to open FITS file pointer or create new file -% (which then contains only the current frequency -% interval) +% outfile - String_Type. If set, write various timing products +% into a FITS file. % elo, ehi - Double_Type: Energy grid, only used for writing outfile % % (Furthermore, the qualifiers of segment_lc_for_psd can be specified) @@ -141,14 +196,10 @@ define lag_energy_spectrum(lc_list, dimseg) { % the lightcurves in lc_list)! The elo, ehi qualifiers are only % for the records in the output file % -% If outfile is a string, a new FITS file will be created, but you -% can also give a file pointer and the function will only create a -% new extension called BAND__Hz. This way you can -% store frequency-resolved lag-energy spectra in one output file by -% always passing the same file pointer (see example). If you want -% to add energy information to the output FITS file, you can pass -% two Double_Type arrays (same length as lc_list) via the elo, ehi -% qualifiers. The columns of the created FITS file are: +% If you want to add energy information to the output FITS file, +% you can pass two Double_Type arrays (same length as lc_list) via +% the elo, ehi qualifiers. The columns of the created extension, +% called BAND__Hz, are: % % ENERGY_LO - Energy channel (keV), taken from elo qualifier % ENERGY_HI - Energy channel (keV), taken from ehi qualifier @@ -175,72 +226,64 @@ define lag_energy_spectrum(lc_list, dimseg) { % channel-of-interest lightcurve % RMS_REF - RMS variability of the reference lightcurve % -% The last three fields should be the same for all lightcurves. -% % If you find any bugs or have questions, please contact ole.koenig@fau.de %\example % %% Example for frequency-resolved lag-energy spectrum -% (flo, fhi) = log_grid(1, 50, n_freqs); % (Hz) +% (f_lo, f_hi) = log_grid(1, 50, n_freqs); % (Hz) % % Energy grid must be the same as extracted lightcurves % (elo, ehi) = log_grid(0.5, 10, n_energies); % (keV) -% variable fp = fits_open_file("test.fits", "c"); -% _for ii (0, n_freqs-1, 1) { % going through frequency range -% variable les = lag_energy_spectrum(lc_names, dimseg -% ; f_min = flo[ii], f_max=fhi[ii], -% verbose=2, outfile=fp, elo=elo, ehi=ehi); -% hplot_with_err(elo, ehi, les.lag, les.errlag); +% variable les = lag_energy_spectrum(lc_names, dimseg +% ; f_lo = f_lo, f_hi=f_hi, +% verbose=2, outfile="test.fits", elo=elo, ehi=ehi); +% _for ii (0, n_freqs-1, 1) { +% ohplot_with_err(elo, ehi, les.lag[*,ii], les.errlag[*,ii]); % } -% fits_close_file(fp); %\seealso{segment_lc_for_psd, foucalc, colacal} %!%- variable dt = qualifier("dt", fits_read_key(lc_list[0], "TIMEDEL")); - variable f_min = qualifier("f_min", 1/(dt*dimseg)); - variable f_max = qualifier("f_max", 1/(2*dt)); % Nyquist frequency + variable f_lo = qualifier("f_lo", [1/(dt*dimseg)]); + variable f_hi = qualifier("f_hi", [1/(2*dt)]); % Nyquist frequency variable verbose = qualifier("verbose", 0); variable deadtime = qualifier("deadtime", 0.0); variable normtype = "Miyamoto"; if (dt == NULL) return "ERROR: Couldn't receive the time resolution from TIMEDEL, need dt qualifier!"; - if (verbose > 0) vmessage("Calculating lag-energy spectrum in frequency range %g-%g Hz", f_min, f_max); variable N_energies = length(lc_list); - variable rmss = Double_Type[N_energies]; - variable rms_refs = Double_Type[N_energies]; - variable lags = Double_Type[N_energies]; - variable errlag_nowak99 = Double_Type[N_energies]; - variable errlag_uttley14 = Double_Type[N_energies]; + variable N_frequencies = length(f_lo); + variable dat = struct{ + lag = Double_Type[N_energies, N_frequencies], + errlag_nowak99 = Double_Type[N_energies, N_frequencies], + errlag_uttley14 = Double_Type[N_energies, N_frequencies], + rms = Double_Type[N_energies], + rms_refs = Double_Type[N_energies], + min_freqs = Double_Type[N_energies, N_frequencies], + max_freqs = Double_Type[N_energies, N_frequencies], + mean_realcpds = Double_Type[N_energies, N_frequencies], + mean_realcpd_errs = Double_Type[N_energies, N_frequencies], + mean_imagcpds = Double_Type[N_energies, N_frequencies], + mean_imagcpd_errs = Double_Type[N_energies, N_frequencies], + mean_phases = Double_Type[N_energies, N_frequencies], + errlag_properr = Double_Type[N_energies, N_frequencies], + errlag_ingram19 = Double_Type[N_energies, N_frequencies], + n_segments = Integer_Type[N_energies, N_frequencies], + n_frequencies = Integer_Type[N_energies, N_frequencies] + }; + variable lc_ref; if (qualifier_exists("lc_ref") == 1){ - lc_ref = getLightcurve(qualifier("lc_ref"), dt, dimseg ;; __qualifiers); + lc_ref = getLightcurve(qualifier("lc_ref"), dt, dimseg;; __qualifiers); } else { lc_ref = lc_list; } - if (qualifier_exists("outfile") == 1){ - variable fp; - if (typeof(qualifier("outfile")) == String_Type){ - fp = fits_open_file(qualifier("outfile"), "c"); - } else { - fp = qualifier("outfile"); - } - variable min_freqs = Double_Type[N_energies]; - variable max_freqs = Double_Type[N_energies]; - variable mean_realcpds = Double_Type[N_energies]; - variable mean_realcpd_errs = Double_Type[N_energies]; - variable mean_imagcpds = Double_Type[N_energies]; - variable mean_imagcpd_errs = Double_Type[N_energies]; - variable mean_phases = Double_Type[N_energies]; - variable errlag_properr = Double_Type[N_energies]; - variable errlag_ingram19 = Double_Type[N_energies]; - variable n_segments = Integer_Type[N_energies]; - variable n_frequencies = Integer_Type[N_energies]; - } - %% Go through the energy-resolved lightcurves (channel-of-interests) - variable ii; - _for ii (0, N_energies-1, 1){ - if (verbose > 1) vmessage(" ***(%s): Calculate time lags for %s", _function_name(), lc_list[ii]); + %% Loop through energy-resolved lightcurves (channel-of-interests) + variable ii, jj; + _for ii (0, N_energies-1, 1){ + + if (verbose > 0) vmessage(" ***(%s): Calculate time lags for %s", _function_name(), lc_list[ii]); variable channel_of_interest = getLightcurve (lc_list[ii], dt, dimseg ;; __qualifiers); variable reference = getReferenceLightcurve (lc_ref, channel_of_interest, ii, dt, dimseg ;; __qualifiers); @@ -250,83 +293,53 @@ define lag_energy_spectrum(lc_list, dimseg) { rate1 = reference.rate, rate2 = channel_of_interest.rate}, dimseg ; normtype=normtype, RMS=&rms, deadtime=deadtime); - (rmss[ii], rms_refs[ii]) = rms[0], rms[1]; - - variable res_f = filterFouQuantForFrequencyRange(res, f_min, f_max ; verbose=verbose); + (dat.rms[ii], dat.rms_refs[ii]) = rms[0], rms[1]; - %% Coherence and lag calculation on averaged cross-spectrum - variable cola = colacal(0.5*(f_min+f_max), % mid frequency - mean(res_f.realcpd12)+mean(res_f.imagcpd12)*1i, - mean(res_f.noicpd12), % not needed for lag+error - mean(res_f.rawpsd1), mean(res_f.rawpsd2), - mean(res_f.noipsd1), mean(res_f.noipsd2), - mean(res_f.sigpsd1), mean(res_f.sigpsd2), % not needed for lag+err - res_f.numavgall[0]*length(res_f.freq)); % K*M with M=number of frequencies averaged over and K=number of segments - (lags[ii], errlag_nowak99[ii], errlag_uttley14[ii]) = cola.lag, cola.errlag, cola.errlag_uttley14; + %% Loop through frequency ranges + _for jj (0, N_frequencies-1, 1){ + if (verbose > 1) vmessage(" ***(%s): Calculating lag-energy spectrum in frequency range %g-%g Hz", + _function_name(), f_lo[jj], f_hi[jj]); + + variable res_f = filterFouQuantForFrequencyRange(res, f_lo[jj], f_hi[jj] ; verbose=verbose); - if (qualifier_exists("outfile") == 1){ + %% Coherence and lag calculation on averaged cross-spectrum + variable cola = colacal(0.5*(f_lo[jj]+f_hi[jj]), % mid frequency + mean(res_f.realcpd12)+mean(res_f.imagcpd12)*1i, + mean(res_f.noicpd12), % not needed for lag+error + mean(res_f.rawpsd1), mean(res_f.rawpsd2), + mean(res_f.noipsd1), mean(res_f.noipsd2), + mean(res_f.sigpsd1), mean(res_f.sigpsd2), % not needed for lag+err + res_f.numavgall[0]*length(res_f.freq)); % K*M with M=number of frequencies averaged over and K=number of segments + (dat.lag[ii,jj], dat.errlag_nowak99[ii,jj], dat.errlag_uttley14[ii,jj]) = + cola.lag, cola.errlag, cola.errlag_uttley14; + variable m_ReCPD = mean(res_f.realcpd12); variable m_ReCPD_err = mean(res_f.errrealcpd12); variable m_ImCPD = mean(res_f.imagcpd12); variable m_ImCPD_err = mean(res_f.errimagcpd12); - + %% Calculation from Gaussian error propagation variable m_phase = atan2(m_ImCPD, m_ReCPD); % radians - errlag_properr[ii]=1/(PI*(f_min+f_max))* + dat.errlag_properr[ii,jj]=1/(PI*(f_lo[jj]+f_hi[jj]))* sqrt((m_ReCPD/(m_ImCPD^2+m_ReCPD^2))^2*m_ImCPD_err^2 + - (m_ImCPD/(m_ImCPD^2+m_ReCPD^2))^2*m_ReCPD_err^2 - ); + (m_ImCPD/(m_ImCPD^2+m_ReCPD^2))^2*m_ReCPD_err^2); - %% Implement error calculation from Ingram 2019, Eq. 19 + %% Error calculation from Ingram 2019, Eq. 19: To be implemented - min_freqs[ii] = min(res_f.freq); - max_freqs[ii] = max(res_f.freq); - mean_realcpds[ii] = m_ReCPD; - mean_realcpd_errs[ii] = m_ReCPD_err; - mean_imagcpds[ii] = m_ImCPD; - mean_imagcpd_errs[ii] = m_ImCPD_err; - mean_phases[ii] = m_phase; - n_segments[ii] = res_f.numavgall[0]; - n_frequencies[ii] = length(res_f.freq); + dat.min_freqs[ii,jj] = min(res_f.freq); + dat.max_freqs[ii,jj] = max(res_f.freq); + dat.mean_realcpds[ii,jj] = m_ReCPD; + dat.mean_realcpd_errs[ii,jj] = m_ReCPD_err; + dat.mean_imagcpds[ii,jj] = m_ImCPD; + dat.mean_imagcpd_errs[ii,jj] = m_ImCPD_err; + dat.mean_phases[ii,jj] = m_phase; + dat.n_segments[ii,jj] = res_f.numavgall[0]; + dat.n_frequencies[ii,jj] = length(res_f.freq); } } - - if (qualifier_exists("outfile") == 1){ - % write_cpd_to_ascii(mean_realcpds, mean_realcpd_errs, mean_imagcpds, mean_imagcpd_errs, f_min, f_max); - - variable extname = sprintf("BAND_%f_%fHz", f_min, f_max); - variable sdata = struct{ - "ENERGY_LO"=qualifier("elo", Double_Type[N_energies]), - "ENERGY_HI"=qualifier("ehi", Double_Type[N_energies]), - "LAG"=lags, - "LAG_ERR_UTTLEY14"=errlag_uttley14, - "LAG_ERR_NOWAK99"=errlag_nowak99, - "LAG_ERR_PROPERR"=errlag_properr, - "LAG_ERR_INGRAM19"=errlag_ingram19, - "MEAN_IMAG_CPD"=mean_imagcpds, - "MEAN_IMAG_CPD_ERR"=mean_imagcpd_errs, - "MEAN_REAL_CPD"=mean_realcpds, - "MEAN_REAL_CPD_ERR"=mean_realcpd_errs, - "MEAN_PHASE"=mean_phases, - "N_SEGMENTS"=n_segments, - "N_FREQUENCIES"=n_frequencies, - "MIN_FREQUENCY"=min_freqs, - "MAX_FREQUENCY"=max_freqs, - "RMS"=rmss, - "RMS_REF"=rms_refs - }; - variable keys = struct{ - F_MIN=f_min, F_MAX=f_max, DIMSEG=dimseg, TIMEDEL=dt, NORMTYPE=normtype - }; - variable hist = struct{ - history = [sprintf("Wrote lag-energy spectrum in frequency range %f-%fHz", f_min, f_max)], - comment = ["Time unit is seconds, energy unit is keV, frequency unit is Hertz", - "Analyzed lightcurves:", lc_list] - }; - fits_write_binary_table(fp, extname, sdata, keys, hist); - - if (typeof(qualifier("outfile")) == String_Type) fits_close_file(fp); - } - return struct{ lag = lags, errlag = errlag_nowak99 }; + if (qualifier_exists("outfile")) write_to_outfile(dat, + f_lo, f_hi, dimseg, dt, normtype, lc_list;; __qualifiers()); + + return dat; } -- GitLab From e837f3edd1b1f2aee2ee5054c8b6de33a89d91a1 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Sat, 5 Feb 2022 01:51:19 +0100 Subject: [PATCH 2/8] foucalc: Add rms qualifier with error calculation - The "RMS" qualifier calculates the RMS of the lightcurves (given by rate1 and rate2), and returns it as an array, however, without error. - A new qualifier "rms" is added, which calculates the RMS of the signal PSD (same as "RMS" qualifier), noise PSD and the error according to Uttley et al., A&A Review 22, 72, 2014, Sect. 2.2.3. - The reference variable is henceforth of Struct_Type - Backwards-compatibility: The "RMS" qualifier is labelled "deprecated" and a warning is is thrown - This commit also remove some old commented-out code --- src/data/timing/foucalc.sl | 54 +++++++++++++++++++++++++++----------- 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/src/data/timing/foucalc.sl b/src/data/timing/foucalc.sl index 1613ec0a..0b2ee941 100644 --- a/src/data/timing/foucalc.sl +++ b/src/data/timing/foucalc.sl @@ -14,8 +14,10 @@ define foucalc() %\qualifier{deadtime}{ [\code{=1e-5}]: detector deadtime in seconds} %\qualifier{fmin}{: minimum frequency used for RMS calculation} %\qualifier{fmax}{: maximum frequency used for RMS calculation} -%\qualifier{RMS}{: reference to a variable to store the RMS of each light curve -% in the [fmin, fmax] band. Can only be used with normtype="Miyamoto".} +%\qualifier{RMS}{: DEPRECATED, use rms qualifier instead.} +%\qualifier{rms}{: reference to a variable to store the signal, noise RMS and error +% of each light curve in the [fmin, fmax] band. Can only be used with +% normtype="Miyamoto".} %\qualifier{avgrate}{: reference to a variable to store the average rate of each light curve} %\qualifier{compact}{: compact the output structure, only keep the most important quantities} %\qualifier{noCPD}{: do not calculate cross power densities and derived quantities} @@ -80,16 +82,22 @@ define foucalc() % Lightcurve parameters variable avgrate = Double_Type[numchan]; - variable rms = Double_Type[numchan]; - + variable rms_deprecated = Double_Type[numchan]; % used for deprecated RMS qualifier + variable fieldnames = ["freq", "numavgall"]; % field names of output structure + variable rms_fieldnames = String_Type[0]; + variable i, j; - _for i (1, numchan, 1) + _for i (1, numchan, 1){ fieldnames = [fieldnames, ["rawpsd", "errpsd", "normpsd", "errnormpsd", "noipsd", "sigpsd", "noinormpsd", "signormpsd", "effnoinormpsd"]+string(i)]; + rms_fieldnames = [rms_fieldnames, ["sigrms", "errrms", "noirms"]+string(i)]; + } + _for i (1, numchan, 1) _for j (i+1, numchan, 1) fieldnames = [fieldnames, ["realcpd", "imagcpd", "errrealcpd", "errimagcpd", "noicpd", "rawcof", "cof", "errcof", "lag", "errlag"]+string(i)+string(j)]; % assuming i < j < 100 ... variable fouquant = @Struct_Type(fieldnames); % create structure + variable rmsquant = @Struct_Type(rms_fieldnames); % Calculating Fourier Frequency Array fouquant.numavgall = ones(dimseg/2) * numseg; @@ -114,7 +122,7 @@ define foucalc() vmessage("warning (foucalc): Deadtime keyword not set. Using default for RXTE: 1e-05 s"); } - % RMS Qualifiers + % RMS Qualifiers (see, e.g., Uttley et al., A&A Review 22, 72, 2014, Sect. 2.2.3) variable fmin = qualifier("fmin", min(fouquant.freq)); variable fmax = qualifier("fmax", max(fouquant.freq)); @@ -157,14 +165,12 @@ define foucalc() % Calculating observational noise with deadtime influence (noipsd, noinormpsd) variable noipsd, noinormpsd; -% (, noipsd) = psdcorr_zhang(avgrate[i-1], tseg, dimseg;; struct_combine(__qualifiers(), struct { avgbkg=avgbkg[i-1], unnormalized })); -% (, noinormpsd) = psdcorr_zhang(avgrate[i-1], tseg, dimseg;; __qualifiers()); (, noipsd, noinormpsd) = psdcorr_zhang(avgrate[i-1], tseg, dimseg;; struct_combine(__qualifiers(), struct { avgbkg=avgbkg[i-1] })); set_struct_field(fouquant, "noipsd"+istr, noipsd); set_struct_field(fouquant, "noinormpsd"+istr, noinormpsd); % effective noise level - set_struct_field(fouquant, "effnoinormpsd"+istr, noinormpsd/sqrt(numseg)); % fouquant.numavgall + set_struct_field(fouquant, "effnoinormpsd"+istr, noinormpsd/sqrt(numseg)); % calculate normalized signal psd variable signormpsd = tempnormpsd - noinormpsd; set_struct_field(fouquant, "signormpsd"+istr, signormpsd); @@ -175,8 +181,15 @@ define foucalc() variable idxmin = min(idx); variable idxmax = max(idx); variable df = fouquant.freq[[idxmin+1:idxmax]] - fouquant.freq[[idxmin:idxmax-1]]; - rms[i-1] = sqrt( sum(signormpsd[[idxmin:idxmax-1]] * df) ); - if(verbose) vmessage("RMS(%g-%g Hz) = %.1f%%", fmin, fmax, rms[i-1]*100); + rms_deprecated[i-1] = sqrt( sum(signormpsd[[idxmin:idxmax-1]] * df) ); + + variable sigrms = sqrt( sum(signormpsd[[idxmin:idxmax-1]] * df) ); + variable noirms = sqrt( sum(noinormpsd[[idxmin:idxmax-1]] * df) ); + variable errrms = sqrt((2.*sigrms^2*noirms^2+noirms^4)/(2.*numseg*length(idx)*sigrms^2)); + set_struct_field(rmsquant, "sigrms"+istr, sigrms); + set_struct_field(rmsquant, "errrms"+istr, errrms); + set_struct_field(rmsquant, "noirms"+istr, noirms); + if(verbose) vmessage("RMS(%g-%g Hz) = (%g +/- %g) %%", fmin, fmax, sigrms*100, errrms*100); } else if(verbose) message(""); % to finish the earlier printf without "\n" @@ -227,19 +240,30 @@ define foucalc() set_struct_field(fouquant, "errcof"+lostr+histr, cola.errcof); set_struct_field(fouquant, "lag" +lostr+histr, cola.lag); set_struct_field(fouquant, "errlag"+lostr+histr, cola.errlag); - %%% set_struct_field(fouquant, "sigcpd"+lostr+histr, sigcpd); cola= NULL; cpd = NULL; errcpd=NULL; } } + variable rmsref = qualifier("rms"); + if(normtype=="Miyamoto" and rmsref!=NULL) + { + if(typeof(rmsref) == Ref_Type) + (@rmsref) = rmsquant; + else + vmessage("warning (%s): rms qualifier has to be a reference", _function_name()); + } + + % Deprecated RMS qualifier (kept for backwards-compatibility) variable RMSref = qualifier("RMS"); if(normtype=="Miyamoto" and RMSref!=NULL) { - if(typeof(RMSref)==Ref_Type) - (@RMSref) = rms; - else + if(typeof(RMSref)==Ref_Type){ + vmessage("warning (%s): RMS qualifier is deprecated, please change to rms.", _function_name()); + (@RMSref) = rms_deprecated; + } else { vmessage("warning (%s): RMS qualifier has to be a reference", _function_name()); + } } variable avgrateref = qualifier("avgrate"); -- GitLab From 09d144056acda8a81ad9f4783bc7ac892c762997 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Sat, 5 Feb 2022 01:57:10 +0100 Subject: [PATCH 3/8] lag_energy_spectrum: Add freq.-resolved RMS calculation - The RMS in each frequency band for each energy bin is now calculated and written into the output file - Calculation is the same as in foucalc. However, there might be a bug in foucalc, namely the RMS frequency filtering in foucalc selects the range [fmin,fmax] but it should be ]fmin,fmax] as otherwise lower bins are counted twice when looped over multiple frequencies. In the lag_energy_spectrum function, the filtering is done on ]fmin,fmax]. This likely has no influence on the RMS calculation as the DFT grid is most likely different from the user's frequency grid. - Write n_segments, n_frequencies, f_min_dft, f_max_dft as keywords instead of columns because they must be the same for each given energy-resolved lightcurve (add check which throws a warning if this is not the case) --- src/data/timing/lag_energy_spectrum.sl | 238 +++++++++++++++---------- 1 file changed, 139 insertions(+), 99 deletions(-) diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl index cdf35b40..4eb9c0a9 100644 --- a/src/data/timing/lag_energy_spectrum.sl +++ b/src/data/timing/lag_energy_spectrum.sl @@ -1,15 +1,20 @@ % -*- mode: slang; mode: fold -*- -private define filterFouQuantForFrequencyRange(fouquan, f_min, f_max) { +private define filterFouQuantForFrequencyRange(fouquant, fmin, fmax) { %{{{ - variable filt = where(f_min < fouquan.freq <= f_max); - if (length(filt) == 0) vmessage("ERROR (lag_energy_spectrum): No frequency bins between %g-%g in Fourier products!", f_min, f_max); + variable idx = where(fmin < fouquant.freq <= fmax); + variable idxmin = min(idx); + variable idxmax = max(idx); + variable df = fouquant.freq[[idxmin+1:idxmax]] - fouquant.freq[[idxmin:idxmax-1]]; + + if (length(idx) == 0) vmessage("ERROR (lag_energy_spectrum): No frequency bins between %g-%g in Fourier products!", fmin, fmax); if (qualifier("verbose", 0) > 1){ - vmessage(" ***(lag_energy_spectrum): Minimal frequency bin %d at %g Hz", filt[0], fouquan.freq[filt[0]]); - vmessage(" ***(lag_energy_spectrum): Maximal frequency bin %d at %g Hz", filt[-1], fouquan.freq[filt[-1]]); - vmessage(" ***(lag_energy_spectrum): Number of frequency bins: %d", length(filt)); + vmessage(" ***(lag_energy_spectrum): Minimal frequency bin %d at %g Hz", idxmin, fouquant.freq[idxmin]); + vmessage(" ***(lag_energy_spectrum): Maximal frequency bin %d at %g Hz", idxmax, fouquant.freq[idxmax]); + vmessage(" ***(lag_energy_spectrum): Number of frequency bins: %d", length(idx)); } - return struct_filter(fouquan, filt ; copy); + variable ret = struct_filter(fouquant, idx; copy); + return struct_combine(ret, struct{df = df, idxmin = idxmin, idxmax = idxmax, n_freqs = length(idx)}); } %}}} private define verifyLightcurveTime(lc1, lc2){ @@ -76,11 +81,11 @@ private define getReferenceLightcurve(lc_ref, cof_lc, index, dt, dimseg){ return struct{time=lc_ref.time, rate=rate}; } %}}} -private define write_cpd_to_ascii(realcpd, realcpd_err, imagcpd, imagcpd_err, f_min, f_max){ +private define write_cpd_to_ascii(realcpd, realcpd_err, imagcpd, imagcpd_err, fmin, fmax){ %{{{ % ASCII files of the cross-spectrum are needed for fitting with RELTRANS - variable fptr_imag = fopen(sprintf("freq_cross_spec_imag_%f_%fHz.ascii", f_min, f_max), "w"); - variable fptr_real = fopen(sprintf("freq_cross_spec_real_%f_%fHz.ascii", f_min, f_max), "w"); + variable fptr_imag = fopen(sprintf("freq_cross_spec_imag_%f_%fHz.ascii", fmin, fmax), "w"); + variable fptr_real = fopen(sprintf("freq_cross_spec_real_%f_%fHz.ascii", fmin, fmax), "w"); variable ii; _for ii (0, length(realcpd)-1, 1){ () = fputs(sprintf("%i\t%f\t%f\n", ii+1, imagcpd[ii], imagcpd_err[ii]), fptr_imag); @@ -90,19 +95,35 @@ private define write_cpd_to_ascii(realcpd, realcpd_err, imagcpd, imagcpd_err, f_ () = fclose(fptr_real); } %}}} -private define write_to_outfile(dat, f_lo, f_hi, dimseg, dt, normtype, lc_list){ +private define verify_frequency_binning(dat){ + %{{{ + %% Check that the frequency binning of every energy-resolved + %% lightcurve was the same for the individual frequency bins. If + %% this is not the case this means that the energy-resolved + %% lightcurves were extracted with different binning/time + %% resolution. + variable ii; + variable N_energies = array_shape(dat.lag)[0]; + _for ii (1, N_energies-1, 1){ + if (_eqs(dat.n_segments[ii,*], dat.n_segments[0,*]) == 0 or + _eqs(dat.n_frequencies[ii,*], dat.n_frequencies[0,*]) == 0 or + _eqs(dat.f_min_dft[ii,*], dat.f_min_dft[0,*]) == 0 or + _eqs(dat.f_max_dft[ii,*], dat.f_max_dft[0,*]) == 0){ + vmessage("WARNING (%s): The DFT binning of the energy-resolved lightcurves differs! Check the time resolutios of your individual lightcurves!", + _function_name()); + } + } +} %}}} + +private define write_to_outfile(dat, f_lo, f_hi, keywords, lc_list){ %{{{ - % To Do: MIN_FREQS, MAX_FREQS, N_FREQUENCIES, and N_SEGMENTS should - % be written as keywords, not columns, because they're the same - % within each frequency extension. - variable fp = fits_open_file(qualifier("outfile"), "c"); - variable N_energies, N_frequencies; - (N_energies, N_frequencies) = array_shape(dat.lag)[0], array_shape(dat.lag)[1]; + variable N_energies, N_les; + (N_energies, N_les) = array_shape(dat.lag)[0], array_shape(dat.lag)[1]; variable jj; - _for jj (0, N_frequencies-1, 1){ + _for jj (0, N_les-1, 1){ % write_cpd_to_ascii(dat.mean_realcpds[*,jj], dat.mean_realcpd_errs[*,jj], % dat.mean_imagcpds[*,jj], dat.mean_imagcpd_errs[*,jj], % f_lo[jj], f_hi[jj]); @@ -115,22 +136,20 @@ private define write_to_outfile(dat, f_lo, f_hi, dimseg, dt, normtype, lc_list){ "LAG_ERR_UTTLEY14" = dat.errlag_uttley14[*,jj], "LAG_ERR_NOWAK99" = dat.errlag_nowak99[*,jj], "LAG_ERR_PROPERR" = dat.errlag_properr[*,jj], - "LAG_ERR_INGRAM19" = dat.errlag_ingram19[*,jj], + % "LAG_ERR_INGRAM19" = dat.errlag_ingram19[*,jj], "MEAN_IMAG_CPD" = dat.mean_imagcpds[*,jj], "MEAN_IMAG_CPD_ERR" = dat.mean_imagcpd_errs[*,jj], "MEAN_REAL_CPD" = dat.mean_realcpds[*,jj], "MEAN_REAL_CPD_ERR" = dat.mean_realcpd_errs[*,jj], "MEAN_PHASE" = dat.mean_phases[*,jj], - "N_SEGMENTS" = dat.n_segments[*,jj], - "N_FREQUENCIES" = dat.n_frequencies[*,jj], - "MIN_FREQUENCY" = dat.min_freqs[*,jj], - "MAX_FREQUENCY" = dat.max_freqs[*,jj], - "RMS" = dat.rms, - "RMS_REF" = dat.rms_refs - }; - variable keys = struct{ - F_MIN=f_lo[jj], F_MAX=f_hi[jj], DIMSEG=dimseg, TIMEDEL=dt, NORMTYPE=normtype + "RMS" = dat.rms[*,jj], + "RMS_ERR" = dat.rms_err[*,jj], }; + + variable keys = struct_combine(struct{ + F_MIN=f_lo[jj], F_MAX=f_hi[jj], N_SEGMENTS=dat.n_segments[0,jj], + N_FREQUENCIES=dat.n_frequencies[0,jj], F_MIN_DFT=dat.f_min_dft[0,jj], + F_MAX_DFT=dat.f_max_dft[0,jj]}, keywords); variable hist = struct{ history = [sprintf("Wrote lag-energy spectrum in frequency range %f-%fHz", f_lo[jj], f_hi[jj])], comment = ["Time unit is seconds, energy unit is keV, frequency unit is Hertz", @@ -140,8 +159,8 @@ private define write_to_outfile(dat, f_lo, f_hi, dimseg, dt, normtype, lc_list){ } fits_close_file(fp); - () = system(sprintf("fthedit %s\[0] keyword=N_FREQS operation=add value=%i comment=\"Number of frequency bands\"", - qualifier("outfile"), N_frequencies)); + () = system(sprintf("fthedit %s\[0] keyword=N_LES operation=add value=%i comment=\"Number of frequency-resolved lag-energy spectra\"", + qualifier("outfile"), N_les)); vmessage("Wrote lag-energy spectrum data to %s", qualifier("outfile")); } %}}} @@ -168,7 +187,7 @@ define lag_energy_spectrum(lc_list, dimseg) { % deadtime - Detector deadtime in seconds [default: 0.0s] % f_lo - Array_Type. For frequency-resolved lag-energy % spectra. Default: lowest sampled frequency [an array -% containing only f_min=1/(dt*dimseg)]. Unit: Hz +% containing only fmin=1/(dt*dimseg)]. Unit: Hz % f_hi - See f_lo. Default: Nyquist frequency [1/(2*dt)]. Unit: Hz % verbose - Increase verbosity [default=0] % outfile - String_Type. If set, write various timing products @@ -187,7 +206,7 @@ define lag_energy_spectrum(lc_list, dimseg) { % * The average time lag in the frequency interval is calculated by % the normal mean on the imaginary and real part of the % cross-spectrum: -% atan2(mean(Imag(CPD)), mean(Real(CPD)))/(PI*(f_min+f_max)) +% atan2(mean(Imag(CPD)), mean(Real(CPD)))/(PI*(fmin+fmax)) % * Reference lightcurve: For each energy bin, the function takes the % summed lightcurves of all but the current energy bands (always % excluding the lightcurve for the cross-spectrum). @@ -199,32 +218,35 @@ define lag_energy_spectrum(lc_list, dimseg) { % If you want to add energy information to the output FITS file, % you can pass two Double_Type arrays (same length as lc_list) via % the elo, ehi qualifiers. The columns of the created extension, -% called BAND__Hz, are: +% called BAND__Hz, are: % -% ENERGY_LO - Energy channel (keV), taken from elo qualifier -% ENERGY_HI - Energy channel (keV), taken from ehi qualifier -% MEAN_PHASE - Phase (radians), computed as atan2(mean(Im(CPD)), mean(Re(CPD))) -% LAG - Time lag (s), computed as phase/(PI*(f_min+f_max)) -% LAG_ERR_UTTLEY+14 - Error computed from coherence with bias -% Uttley et al., 2014, A&A Rv. 22, 72 (Eq. 11-12) -% LAG_ERR_NOWAK+99 - Error computed from coherence without bias (default) -% Nowak et al., 1999, ApJ, 510, 874 (Sect. 4) -% LAG_ERR_PROPERR - Error computed from Gaussian error propagation -% LAG_ERR_INGRAM19 - (to be implemented) -% MEAN_REAL/IMAG_CPD - real/imaginary part of the mean of the -% segment-averaged cross-spectrum +% ENERGY_LO - Energy channel (keV), taken from elo qualifier +% ENERGY_HI - Energy channel (keV), taken from ehi qualifier +% MEAN_PHASE - Phase (radians), computed as atan2(mean(Im(CPD)), mean(Re(CPD))) +% LAG - Time lag (s), computed as phase/(PI*(fmin+fmax)) +% LAG_ERR_UTTLEY+14 - Error computed from coherence with bias +% Uttley et al., 2014, A&A Rv. 22, 72 (Eq. 11-12) +% LAG_ERR_NOWAK+99 - Error computed from coherence without bias +% Nowak et al., 1999, ApJ, 510, 874 (Sect. 4) +% LAG_ERR_PROPERR - Error computed from Gaussian error propagation +% MEAN_REAL/IMAG_CPD - real/imaginary part of the mean of the +% segment-averaged cross-spectrum % MEAN_REAL/IMAG_CPD_ERR - standard error on the mean of the -% real/imaginary segment-averaged cross-spectrum -% N_SEGMENTS - Number of segments averaged over (foucalc.numavgall) -% N_FREQUENCIES - Number of frequencies averaged over (as in -% "mean of the CPD") -% MIN/MAX_FREQUENCY - Minimal/Maximal frequency of the frequency-filtered -% Fourier products - this can be different from -% f_min due to the discretization of the -% *Discrete* Fourier Transform -% RMS - Root Mean Squared variability of the -% channel-of-interest lightcurve -% RMS_REF - RMS variability of the reference lightcurve +% real/imaginary segment-averaged cross-spectrum +% RMS/RMS_ERR - Root Mean Squared variability of the +% channel-of-interest lightcurve +% +% Following keywords are written into each extension: +% +% N_SEGMENTS - Number of segments averaged over (foucalc.numavgall) +% N_FREQUENCIES - Number of frequencies averaged over (as in +% "mean of the CPD") +% F_MIN/F_MAX - Minimum and maximum frequency as given by +% f_lo, f_hi qualifiers +% F_MIN_DFT/F_MAX_DFT - Minimal/Maximal frequency of the frequency-filtered +% Fourier products - this can be different from +% F_MIN/F_MAX due to the discretization of the +% *Discrete* Fourier Transform % % If you find any bugs or have questions, please contact ole.koenig@fau.de %\example @@ -250,25 +272,27 @@ define lag_energy_spectrum(lc_list, dimseg) { if (dt == NULL) return "ERROR: Couldn't receive the time resolution from TIMEDEL, need dt qualifier!"; variable N_energies = length(lc_list); - variable N_frequencies = length(f_lo); + variable N_les = length(f_lo); variable dat = struct{ - lag = Double_Type[N_energies, N_frequencies], - errlag_nowak99 = Double_Type[N_energies, N_frequencies], - errlag_uttley14 = Double_Type[N_energies, N_frequencies], - rms = Double_Type[N_energies], - rms_refs = Double_Type[N_energies], - min_freqs = Double_Type[N_energies, N_frequencies], - max_freqs = Double_Type[N_energies, N_frequencies], - mean_realcpds = Double_Type[N_energies, N_frequencies], - mean_realcpd_errs = Double_Type[N_energies, N_frequencies], - mean_imagcpds = Double_Type[N_energies, N_frequencies], - mean_imagcpd_errs = Double_Type[N_energies, N_frequencies], - mean_phases = Double_Type[N_energies, N_frequencies], - errlag_properr = Double_Type[N_energies, N_frequencies], - errlag_ingram19 = Double_Type[N_energies, N_frequencies], - n_segments = Integer_Type[N_energies, N_frequencies], - n_frequencies = Integer_Type[N_energies, N_frequencies] + lag = Double_Type[N_energies, N_les], + errlag_nowak99 = Double_Type[N_energies, N_les], + errlag_uttley14 = Double_Type[N_energies, N_les], + rms = Double_Type[N_energies, N_les], + rms_err = Double_Type[N_energies, N_les], + min_freqs = Double_Type[N_energies, N_les], + max_freqs = Double_Type[N_energies, N_les], + mean_realcpds = Double_Type[N_energies, N_les], + mean_realcpd_errs = Double_Type[N_energies, N_les], + mean_imagcpds = Double_Type[N_energies, N_les], + mean_imagcpd_errs = Double_Type[N_energies, N_les], + mean_phases = Double_Type[N_energies, N_les], + errlag_properr = Double_Type[N_energies, N_les], + % errlag_ingram19 = Double_Type[N_energies, N_les], + n_segments = Integer_Type[N_energies, N_les], + n_frequencies = Integer_Type[N_energies, N_les], + f_min_dft = Double_Type[N_energies, N_les], + f_max_dft = Double_Type[N_energies, N_les] }; variable lc_ref; @@ -286,37 +310,49 @@ define lag_energy_spectrum(lc_list, dimseg) { if (verbose > 0) vmessage(" ***(%s): Calculate time lags for %s", _function_name(), lc_list[ii]); variable channel_of_interest = getLightcurve (lc_list[ii], dt, dimseg ;; __qualifiers); - variable reference = getReferenceLightcurve (lc_ref, channel_of_interest, ii, dt, dimseg ;; __qualifiers); + variable reference = getReferenceLightcurve (lc_ref, channel_of_interest, ii, dt, dimseg;; __qualifiers); - variable rms; - variable res = foucalc( struct{ time = channel_of_interest.time, - rate1 = reference.rate, - rate2 = channel_of_interest.rate}, - dimseg ; normtype=normtype, RMS=&rms, deadtime=deadtime); - (dat.rms[ii], dat.rms_refs[ii]) = rms[0], rms[1]; + variable fouquant = foucalc(struct{time = channel_of_interest.time, + rate1 = reference.rate, + rate2 = channel_of_interest.rate}, + dimseg; normtype=normtype, deadtime=deadtime); %% Loop through frequency ranges - _for jj (0, N_frequencies-1, 1){ + _for jj (0, N_les-1, 1){ if (verbose > 1) vmessage(" ***(%s): Calculating lag-energy spectrum in frequency range %g-%g Hz", _function_name(), f_lo[jj], f_hi[jj]); - variable res_f = filterFouQuantForFrequencyRange(res, f_lo[jj], f_hi[jj] ; verbose=verbose); - + variable fouquant_f = filterFouQuantForFrequencyRange(fouquant, f_lo[jj], f_hi[jj]; verbose=verbose); + + % RMS calculation (same as in foucalc but we want it for every + % energy lightcurve and frequency range - running foucalc + % everytime drastically increases the runtime due to the DFT!) + % FIX ME: should only use fouqant_f but somehow df is 1 too short!? + variable sigrms = sqrt( sum(fouquant.signormpsd2[[fouquant_f.idxmin:fouquant_f.idxmax-1]] * fouquant_f.df) ); + variable noirms = sqrt( sum(fouquant.noinormpsd2[[fouquant_f.idxmin:fouquant_f.idxmax-1]] * fouquant_f.df) ); + variable errrms = sqrt((2.*sigrms^2*noirms^2+noirms^4)/(2.*fouquant_f.numavgall[0]*fouquant_f.n_freqs*sigrms^2)); + % variable sigrms = sqrt( sum(fouquant_f.signormpsd2 * fouquant_f.df) ); + % variable noirms = sqrt( sum(fouquant_f.noinormpsd2 * fouquant_f.df) ); + % variable errrms = sqrt((2.*sigrms^2*noirms^2+noirms^4)/(2.*fouquant_f.numavgall[0]*fouquant_f.n_freqs*sigrms^2)); + + dat.rms[ii,jj] = sigrms; + dat.rms_err[ii,jj] = errrms; + %% Coherence and lag calculation on averaged cross-spectrum variable cola = colacal(0.5*(f_lo[jj]+f_hi[jj]), % mid frequency - mean(res_f.realcpd12)+mean(res_f.imagcpd12)*1i, - mean(res_f.noicpd12), % not needed for lag+error - mean(res_f.rawpsd1), mean(res_f.rawpsd2), - mean(res_f.noipsd1), mean(res_f.noipsd2), - mean(res_f.sigpsd1), mean(res_f.sigpsd2), % not needed for lag+err - res_f.numavgall[0]*length(res_f.freq)); % K*M with M=number of frequencies averaged over and K=number of segments + mean(fouquant_f.realcpd12)+mean(fouquant_f.imagcpd12)*1i, + mean(fouquant_f.noicpd12), % not needed for lag+error + mean(fouquant_f.rawpsd1), mean(fouquant_f.rawpsd2), + mean(fouquant_f.noipsd1), mean(fouquant_f.noipsd2), + mean(fouquant_f.sigpsd1), mean(fouquant_f.sigpsd2), % not needed for lag+err + fouquant_f.numavgall[0]*length(fouquant_f.freq)); % K*M with M=number of frequencies averaged over and K=number of segments (dat.lag[ii,jj], dat.errlag_nowak99[ii,jj], dat.errlag_uttley14[ii,jj]) = cola.lag, cola.errlag, cola.errlag_uttley14; - variable m_ReCPD = mean(res_f.realcpd12); - variable m_ReCPD_err = mean(res_f.errrealcpd12); - variable m_ImCPD = mean(res_f.imagcpd12); - variable m_ImCPD_err = mean(res_f.errimagcpd12); + variable m_ReCPD = mean(fouquant_f.realcpd12); + variable m_ReCPD_err = mean(fouquant_f.errrealcpd12); + variable m_ImCPD = mean(fouquant_f.imagcpd12); + variable m_ImCPD_err = mean(fouquant_f.errimagcpd12); %% Calculation from Gaussian error propagation variable m_phase = atan2(m_ImCPD, m_ReCPD); % radians @@ -326,20 +362,24 @@ define lag_energy_spectrum(lc_list, dimseg) { %% Error calculation from Ingram 2019, Eq. 19: To be implemented - dat.min_freqs[ii,jj] = min(res_f.freq); - dat.max_freqs[ii,jj] = max(res_f.freq); dat.mean_realcpds[ii,jj] = m_ReCPD; dat.mean_realcpd_errs[ii,jj] = m_ReCPD_err; dat.mean_imagcpds[ii,jj] = m_ImCPD; dat.mean_imagcpd_errs[ii,jj] = m_ImCPD_err; dat.mean_phases[ii,jj] = m_phase; - dat.n_segments[ii,jj] = res_f.numavgall[0]; - dat.n_frequencies[ii,jj] = length(res_f.freq); + dat.f_min_dft[ii,jj] = min(fouquant_f.freq); + dat.f_max_dft[ii,jj] = max(fouquant_f.freq); + dat.n_segments[ii,jj] = fouquant_f.numavgall[0]; + dat.n_frequencies[ii,jj] = length(fouquant_f.freq); } } + + verify_frequency_binning(dat); - if (qualifier_exists("outfile")) write_to_outfile(dat, - f_lo, f_hi, dimseg, dt, normtype, lc_list;; __qualifiers()); + if (qualifier_exists("outfile")){ + variable keywords = struct{DIMSEG=dimseg, TIMEDEL=dt, NORMTYPE=normtype}; + write_to_outfile(dat, f_lo, f_hi, keywords, lc_list;; __qualifiers()); + } return dat; } -- GitLab From 55d436ac4270f7a27d917798ac43234011252674 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Fri, 18 Feb 2022 19:48:47 +0100 Subject: [PATCH 4/8] lag_energy_spectrum: Write PSD of each energy channel into outfile --- src/data/timing/lag_energy_spectrum.sl | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl index 4eb9c0a9..61f99474 100644 --- a/src/data/timing/lag_energy_spectrum.sl +++ b/src/data/timing/lag_energy_spectrum.sl @@ -121,7 +121,7 @@ private define write_to_outfile(dat, f_lo, f_hi, keywords, lc_list){ variable N_energies, N_les; (N_energies, N_les) = array_shape(dat.lag)[0], array_shape(dat.lag)[1]; - + variable jj; _for jj (0, N_les-1, 1){ % write_cpd_to_ascii(dat.mean_realcpds[*,jj], dat.mean_realcpd_errs[*,jj], @@ -158,6 +158,16 @@ private define write_to_outfile(dat, f_lo, f_hi, keywords, lc_list){ fits_write_binary_table(fp, extname, sdata, keys, hist); } + %% Write power spectrum into "PSD" extension + variable colnames = String_Type[N_energies+1]; + colnames[0] = "freq"; + variable ii; + _for ii (0, N_energies-1, 1) colnames[ii+1] = sprintf("psd%i", ii); + variable psddata = @Struct_Type(colnames); + set_struct_field(psddata, "freq", dat.freq[0]); % frequency grids are all the same (because LCs have same binning) + _for ii (0, N_energies-1, 1) set_struct_field(psddata, colnames[ii+1], dat.psd[ii]); + fits_write_binary_table(fp, "PSD", psddata); + fits_close_file(fp); () = system(sprintf("fthedit %s\[0] keyword=N_LES operation=add value=%i comment=\"Number of frequency-resolved lag-energy spectra\"", qualifier("outfile"), N_les)); @@ -288,7 +298,9 @@ define lag_energy_spectrum(lc_list, dimseg) { mean_imagcpd_errs = Double_Type[N_energies, N_les], mean_phases = Double_Type[N_energies, N_les], errlag_properr = Double_Type[N_energies, N_les], - % errlag_ingram19 = Double_Type[N_energies, N_les], + % errlag_ingram19 = Double_Type[N_energies, N_les], + freq = Array_Type[N_energies], + psd = Array_Type[N_energies], n_segments = Integer_Type[N_energies, N_les], n_frequencies = Integer_Type[N_energies, N_les], f_min_dft = Double_Type[N_energies, N_les], @@ -302,7 +314,6 @@ define lag_energy_spectrum(lc_list, dimseg) { lc_ref = lc_list; } - %% Loop through energy-resolved lightcurves (channel-of-interests) variable ii, jj; _for ii (0, N_energies-1, 1){ @@ -315,7 +326,9 @@ define lag_energy_spectrum(lc_list, dimseg) { variable fouquant = foucalc(struct{time = channel_of_interest.time, rate1 = reference.rate, rate2 = channel_of_interest.rate}, - dimseg; normtype=normtype, deadtime=deadtime); + dimseg; normtype=normtype, deadtime=deadtime); + dat.freq[ii] = fouquant.freq; + dat.psd[ii] = fouquant.signormpsd2; %% Loop through frequency ranges _for jj (0, N_les-1, 1){ -- GitLab From e6e1c6662d2df17bbd5942e4e6773fcf4bc30315 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Sat, 19 Feb 2022 00:15:01 +0100 Subject: [PATCH 5/8] foucalc: Fix bug in usage of rate in input struct * Previously, foucalc required the strict format struct{time, rate1[, rate2] as input. This was rather annoying, because often lightcurves have the colums {time, rate, error}. * Commit 56dfa43d implemented the usage of rate instead of rate1. However, this worked only for a struct with the fields {time, rate} or - if more fields were used - one had to set the noCPD qualifier * This is related to the way how the foucalc code infers the numbers of lightcurves it should process. Because of the rigit structure, the code used the length of the struct minus 1. This crashes as soon as there are additional fields in the struct * This commit only counts the number of fields with "rate" in their name to deduce the number of lightcurves to be processed * Also disable support of "rate" usage for CPD calculation: it can be confusing, and the user should be required to use rate1 and rate2 --- src/data/timing/foucalc.sl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/src/data/timing/foucalc.sl b/src/data/timing/foucalc.sl index 0b2ee941..a1a4588e 100644 --- a/src/data/timing/foucalc.sl +++ b/src/data/timing/foucalc.sl @@ -28,7 +28,8 @@ define foucalc() % \code{lc = struct { time=[t_1, t_2, ...], rate1=[r1_1, r1_2, ...], rate2=... };}. % [However, one can also use an array of structures (with an enormous overhead), % \code{lc = [ struct { time=t_1, rate1=r1_1, rate2=...}, struct { time=t_2, rate1=r1_2, rate2=... }, ... ];}.] -% Also specifying "rate" instead of "rate1" is possible. +% Also specifying "rate" instead of "rate1" is possible if no CPD +% should be computed. % % \code{dimseg} is the segment size used for the FFTs, which should therefore be a power of 2. % @@ -78,7 +79,8 @@ define foucalc() variable numseg = int(dimlc/dimseg); variable tseg = bintime * dimseg; variable timelc = bintime * dimlc; - variable numchan = length(get_struct_field_names(lc[0])) - 1; % lc[0] = lc, if lc is no array + variable str_fields = get_struct_field_names(lc[0]); % lc[0]=lc, if lc is no array but struct + variable numchan = length(str_fields[where(is_substr(str_fields, "rate"))]); % Lightcurve parameters variable avgrate = Double_Type[numchan]; @@ -213,11 +215,10 @@ define foucalc() if (verbose) vmessage("Calculating cross power spectrum (=> coherence and timelags) for rate%d and rate%d.", lo, hi); - % If "rate1" is not available in struct, use "rate" - variable rate1 = struct_field_exists(lc, "rate"+lostr) ? array_struct_field(lc, "rate"+lostr) : array_struct_field(lc, "rate"); - variable rate2 = array_struct_field(lc, "rate"+histr); variable cpd, errcpd; - (cpd, errcpd) = cross_power_density(rate1, rate2, numseg, dimseg); + (cpd, errcpd) = cross_power_density(array_struct_field(lc, "rate"+lostr), + array_struct_field(lc, "rate"+histr), + numseg, dimseg); set_struct_field(fouquant, "realcpd"+lostr+histr, Real(cpd)); set_struct_field(fouquant, "imagcpd"+lostr+histr, Imag(cpd)); set_struct_field(fouquant, "errrealcpd"+lostr+histr, Real(errcpd)); -- GitLab From 10c11aa55d9f92035ecb0d6d29b1e51cd4da7ed7 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Fri, 11 Mar 2022 00:20:01 +0100 Subject: [PATCH 6/8] foucalc: Remove Uttley+14 error calculation in colacal.sl * Uttley+14 calculate the coherence with the full power spectrum (noise + signal) in the denominator. This calculation was added in commit 4a64f55c. * This calculation differs from the one in Nowak+99 Eq. 16., Ingram+19 Eq. 21, or Vaughan&Nowak 97 Eq. 8 who divide by the SIGNAL power spectrum only. * Thus, the Uttley+14 coherence will lead to a *smaller* value and thus bigger uncertainties on the time lag. That is, Uttley+14 *overestimate* the uncertainties. * This commit reverses 4a64f55c: It removes the Uttley+14 calculation. --- src/data/timing/colacal.sl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/data/timing/colacal.sl b/src/data/timing/colacal.sl index a865f533..53437a58 100644 --- a/src/data/timing/colacal.sl +++ b/src/data/timing/colacal.sl @@ -23,7 +23,6 @@ define colacal(freq, cpd, noicpd, lopsd, hipsd, noilopsd, noihipsd, siglopsd, si % errcof - one-sigma uncertainty of cof [Vaughan & Nowak, 1997, ApJ, 474, L43 (Eqn. 8, Part 2)]\n % lag - time lag [Nowak et al., 1999, ApJ, 510, 874 (Sect. 4)]\n % errlag - one-sigma uncertainty of lag [Nowak et al., 1999, ApJ, 510, 874 (Eqn. 16)]\n -% errlag_Uttley14 - one-sigma uncertainty of lag [Uttley et al., 2014, A&A Rv. 22, 72 (Eq. 11-12)] % sigcpd - noise-corrected cross-power-density [sigcpd = cpd - noicpd]\n %!%- { @@ -36,9 +35,6 @@ define colacal(freq, cpd, noicpd, lopsd, hipsd, noilopsd, noihipsd, siglopsd, si ); variable lag = atan2(Imag(cpd), Real(cpd)) / (2.*PI*freq); variable errlag = sqrt( (1.-rawcof)/(2.*rawcof)/alln ) / (2.*PI*freq); - variable rawcof_w_bias = (sqr(abs(cpd)) - noicpd) / (lopsd * hipsd); % Uttley+14 Eq. 11, note that this divides by psd and not sigpsd - variable errlag_uttley14 = sqrt( (1.-rawcof_w_bias)/(2.*rawcof)/alln ) / (2.*PI*freq); % Uttley+14 Eq. 12 return struct{rawcof=rawcof, cof=cof, errcof=errcof, - lag=lag, errlag=errlag, errlag_uttley14=errlag_uttley14, - sigcpd=sigcpd}; + lag=lag, errlag=errlag, sigcpd=sigcpd}; } -- GitLab From 87f386070807e4ac6a2d509743bacc886b43d85e Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Sat, 12 Mar 2022 01:48:24 +0100 Subject: [PATCH 7/8] lag_energy_spectrum: Frequency interval filtering from ]fmin,fmax] to [fmin,fmax] * The lower bin should also be included, otherwise the minimal frequency bin (1/T) is never included in the calculation Further notes: * lag_energy_spectrum (and foucalc) filters the frequency range [fmin,fmax] and calculates the RMS in that range with the where function. This can potentially lead to underestimates of the RMS because bins where fmin/fmax is in between bin_lo and bin_hi are NOT counted. * Example: fmin fmax User defined frequency filtering: | | Frequency grid coming from DFT: | | | | | X X only the two bins marked with X are accounted for in the RMS calculation. * A potentially better implementation is to use an if statement to check whether most of the bin is below or above the frequency threshold and then decide to in-/exclude it. --- src/data/timing/lag_energy_spectrum.sl | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl index 61f99474..9e4fa3d5 100644 --- a/src/data/timing/lag_energy_spectrum.sl +++ b/src/data/timing/lag_energy_spectrum.sl @@ -2,7 +2,7 @@ private define filterFouQuantForFrequencyRange(fouquant, fmin, fmax) { %{{{ - variable idx = where(fmin < fouquant.freq <= fmax); + variable idx = where(fmin <= fouquant.freq <= fmax); variable idxmin = min(idx); variable idxmax = max(idx); variable df = fouquant.freq[[idxmin+1:idxmax]] - fouquant.freq[[idxmin:idxmax-1]]; @@ -340,13 +340,9 @@ define lag_energy_spectrum(lc_list, dimseg) { % RMS calculation (same as in foucalc but we want it for every % energy lightcurve and frequency range - running foucalc % everytime drastically increases the runtime due to the DFT!) - % FIX ME: should only use fouqant_f but somehow df is 1 too short!? variable sigrms = sqrt( sum(fouquant.signormpsd2[[fouquant_f.idxmin:fouquant_f.idxmax-1]] * fouquant_f.df) ); variable noirms = sqrt( sum(fouquant.noinormpsd2[[fouquant_f.idxmin:fouquant_f.idxmax-1]] * fouquant_f.df) ); variable errrms = sqrt((2.*sigrms^2*noirms^2+noirms^4)/(2.*fouquant_f.numavgall[0]*fouquant_f.n_freqs*sigrms^2)); - % variable sigrms = sqrt( sum(fouquant_f.signormpsd2 * fouquant_f.df) ); - % variable noirms = sqrt( sum(fouquant_f.noinormpsd2 * fouquant_f.df) ); - % variable errrms = sqrt((2.*sigrms^2*noirms^2+noirms^4)/(2.*fouquant_f.numavgall[0]*fouquant_f.n_freqs*sigrms^2)); dat.rms[ii,jj] = sigrms; dat.rms_err[ii,jj] = errrms; -- GitLab From 4f5fd2644a671c64d9b97fb60854f8ed1d7491d4 Mon Sep 17 00:00:00 2001 From: Ole Koenig Date: Sat, 12 Mar 2022 01:55:54 +0100 Subject: [PATCH 8/8] lag_energy_spectrum: Remove Uttley+14 lag error field * Commit 10c11aa removes the time lag error calculation from Uttley+14 as it only differs from Nowak+99 in the normalization (signal PSD versus full PSD). * Consequently, also the field in the lag-energy spectrum output file needs to be removed --- src/data/timing/lag_energy_spectrum.sl | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/src/data/timing/lag_energy_spectrum.sl b/src/data/timing/lag_energy_spectrum.sl index 9e4fa3d5..851a94a6 100644 --- a/src/data/timing/lag_energy_spectrum.sl +++ b/src/data/timing/lag_energy_spectrum.sl @@ -133,9 +133,8 @@ private define write_to_outfile(dat, f_lo, f_hi, keywords, lc_list){ "ENERGY_LO" = qualifier("elo", Double_Type[N_energies]), "ENERGY_HI" = qualifier("ehi", Double_Type[N_energies]), "LAG" = dat.lag[*,jj], - "LAG_ERR_UTTLEY14" = dat.errlag_uttley14[*,jj], - "LAG_ERR_NOWAK99" = dat.errlag_nowak99[*,jj], - "LAG_ERR_PROPERR" = dat.errlag_properr[*,jj], + "LAG_ERR" = dat.errlag[*,jj], + % "LAG_ERR_PROPERR" = dat.errlag_properr[*,jj], % "LAG_ERR_INGRAM19" = dat.errlag_ingram19[*,jj], "MEAN_IMAG_CPD" = dat.mean_imagcpds[*,jj], "MEAN_IMAG_CPD_ERR" = dat.mean_imagcpd_errs[*,jj], @@ -208,7 +207,7 @@ define lag_energy_spectrum(lc_list, dimseg) { %} %\notes % * The default error computation is Nowak et al., 1999, ApJ, 510, -% 874 (Sect. 4), i.e., from coherence without bias factor. +% 874 (Sect. 4). % * All input lightcurves are accounted for their gaps and segmented % by segment_lc_for_psd. % * The cross-spectrum (CPD) and lags are calculated with foucalc. @@ -234,10 +233,7 @@ define lag_energy_spectrum(lc_list, dimseg) { % ENERGY_HI - Energy channel (keV), taken from ehi qualifier % MEAN_PHASE - Phase (radians), computed as atan2(mean(Im(CPD)), mean(Re(CPD))) % LAG - Time lag (s), computed as phase/(PI*(fmin+fmax)) -% LAG_ERR_UTTLEY+14 - Error computed from coherence with bias -% Uttley et al., 2014, A&A Rv. 22, 72 (Eq. 11-12) -% LAG_ERR_NOWAK+99 - Error computed from coherence without bias -% Nowak et al., 1999, ApJ, 510, 874 (Sect. 4) +% LAG_ERR - Error computed as in Nowak et al., 1999, ApJ, 510, 874 Eq. 16 % LAG_ERR_PROPERR - Error computed from Gaussian error propagation % MEAN_REAL/IMAG_CPD - real/imaginary part of the mean of the % segment-averaged cross-spectrum @@ -286,8 +282,7 @@ define lag_energy_spectrum(lc_list, dimseg) { variable dat = struct{ lag = Double_Type[N_energies, N_les], - errlag_nowak99 = Double_Type[N_energies, N_les], - errlag_uttley14 = Double_Type[N_energies, N_les], + errlag = Double_Type[N_energies, N_les], rms = Double_Type[N_energies, N_les], rms_err = Double_Type[N_energies, N_les], min_freqs = Double_Type[N_energies, N_les], @@ -355,8 +350,7 @@ define lag_energy_spectrum(lc_list, dimseg) { mean(fouquant_f.noipsd1), mean(fouquant_f.noipsd2), mean(fouquant_f.sigpsd1), mean(fouquant_f.sigpsd2), % not needed for lag+err fouquant_f.numavgall[0]*length(fouquant_f.freq)); % K*M with M=number of frequencies averaged over and K=number of segments - (dat.lag[ii,jj], dat.errlag_nowak99[ii,jj], dat.errlag_uttley14[ii,jj]) = - cola.lag, cola.errlag, cola.errlag_uttley14; + (dat.lag[ii,jj], dat.errlag[ii,jj]) = cola.lag, cola.errlag; variable m_ReCPD = mean(fouquant_f.realcpd12); variable m_ReCPD_err = mean(fouquant_f.errrealcpd12); -- GitLab