We can go about doing the task in MATLAB as you described but first, we need to assume some data. Here, I am considering the frequency response of a band pass filter:
%% Sample frequency data
f = linspace(1e9, 2e9, 1000); % Frequency range = 1GHz to 2 GHz
%% Sample data for a band-pass filter
f0 = 1.5e9; % Center frequency at 1.5 GHz
BW = 0.2e9; % Bandwidth of 200 MHz
S21_mag = exp(-((f - f0) / (BW / 2)).^2);
S21_phase = -2 * pi * (f - f0) / max(f);
S21 = S21_mag .* exp(1i * S21_phase);
You may load your own data into variables f and S21 instead of the generic bandpass filter that I have used. The following steps calculate the Correlation function R:
R = @(chi) arrayfun(@(c) trapz(f, S21 .* conj(interp1(f, S21, f + c, 'linear', 0))), chi);
Here R is a function handle that takes an input vector chi and maps each element of chi, say x, to a numeric integral that evaluates R(x). This is done using a function trapz which calculates the numeric integral.
Next, we find the region of chi that correspond to amplitude greater than or equal to ρ|R(0)|
R0 = R(0); % Max value of R
rho = 0.5; % Example level, adjust as needed
chi_values = linspace(-max(f) + min(f), max(f) - min(f), 2000); % Possible χ values (both negative and positive)
R_chi = abs(R(chi_values)); % Calculate R(χ) for each χ
% This finds the chi indices that fall within the bandwidth
B_rho_indices = find(R_chi >= rho * abs(R0));
B_rho = chi_values([B_rho_indices(1), B_rho_indices(end)]); % Corresponding χ values for bandwidth edges
Now, we print the value of the bandwidth along with a plot that makes it easy to visualize:
%% Plotting the outputs:
% Plot the correlation function |R(χ)|
figure;
subplot(2, 1, 1);
plot(chi_values, R_chi, 'b', 'LineWidth', 2);
hold on;
yline(rho * abs(R0), 'r--', 'LineWidth', 2); % Line for ρ|R(0)|
xline(B_rho(1), 'g--', 'LineWidth', 2); % Line for negative Bρ
xline(B_rho(2), 'g--', 'LineWidth', 2); % Line for positive Bρ
xlabel('Frequency Shift χ (Hz)');
ylabel('|R(χ)|');
title('Correlation Function vs Frequency Shift');
legend('|R(χ)|', 'ρ|R(0)|', 'Bρ');
grid on;
% Display the coherence bandwidth
fprintf('Coherence bandwidth Bρ at level ρ = %.2f is from %.2f Hz to %.2f Hz\n', rho, B_rho(1), B_rho(2));
fprintf('Total coherence bandwidth Bρ is %.2f Hz\n', B_rho(2) - B_rho(1));
% Plot the transfer function of the band-pass filter
subplot(2, 1, 2);
plot(f, 20*log10(abs(S21)), 'b', 'LineWidth', 2);
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Transfer Function of the Band-Pass Filter');
grid on;
Along with the output we also get a plot showing the bandwidth between the green dashed lines.
Refer to the following documentation for more information:
Function handles:
trapz:
arrayfun:
