% Generate a signal with white noise and a sinusoidal component
t = linspace(0, 60, 1000);  % Assuming 1000 samples over 60 seconds
white_noise = randn(size(t));
frequency = 0.1;  % Adjust this frequency based on your needs
sinusoidal_component = sin(2*pi*frequency*t);
signal = white_noise + sinusoidal_component;

% Perform Rainflow counting
cycles = rainflow(signal);

% Display the number of cycles
num_cycles = size(cycles, 1);
disp(['Estimated Number of Cycles (Rainflow Counting): ' num2str(num_cycles)]);

% Plot the original signal
figure;
subplot(3, 1, 1);
plot(t, signal);
title('Original Signal');
xlabel('Time (seconds)');
ylabel('Amplitude');

% Plot the identified peaks
subplot(3, 1, 2);
threshold = 0.5;
[peaks, peak_locations] = findpeaks(signal, 'MinPeakHeight', threshold);
plot(t, signal);
hold on;
plot(t(peak_locations), peaks, 'ro', 'MarkerSize', 8);
title('Identified Peaks');
xlabel('Time (seconds)');
ylabel('Amplitude');

% Plot the signal with marked cycles
subplot(3, 1, 3);
plot(t, signal);
hold on;
for i = 1:num_cycles
    start_index = cycles(i, 1);
    end_index = cycles(i, 2);
    
    % Ensure indices are within bounds
    start_index = max(1, start_index);
    end_index = min(length(t), end_index);
    
    plot(t(start_index:end_index), signal(start_index:end_index), 'ro-', 'LineWidth', 2);
end
title('Signal with Rainflow Counting Cycles');
xlabel('Time (seconds)');
ylabel('Amplitude');

% Adjust subplot spacing
subplot(3, 1, 1);
ylim([-5 5]);
subplot(3, 1, 2);
ylim([-5 5]);
subplot(3, 1, 3);
ylim([-5 5]);

hold off;
% ...

% Linear interpolation for smoother surface
[xq, yq] = meshgrid(linspace(min(t), max(t), 500), linspace(min(signal), max(signal), 500));
vq = griddata(time_matrix(:), amplitude_matrix(:), cycle_count_matrix(:), xq, yq);

% Create the 3D surface plot
surf(xq, yq, vq, 'EdgeColor', 'none');
title('Rainflow Counting - 3D Surface Plot');
xlabel('Time (seconds)');
ylabel('Amplitude');
zlabel('Cycle Count');
view(3); % Set the view to 3D

% Additional settings for better visualization
colormap(jet);
colorbar;
hold on;

% Plot the identified cycles on top of the surface plot
for i = 1:num_cycles
    start_index = round(cycles(i, 1));
    end_index = round(cycles(i, 2));
    
    % Ensure indices are within bounds
    start_index = max(1, start_index);
    end_index = min(length(t), end_index);
    
    scatter3(t(start_index:end_index), signal(start_index:end_index), i * ones(size(start_index:end_index)), 'filled', 'r');
end

hold off;

% ...
%%
% Plot the original signal
figure;
subplot(4, 1, 1);
plot(t, signal);
title('Original Signal');
xlabel('Time (seconds)');
ylabel('Amplitude');

% Plot the identified peaks
subplot(4, 1, 2);
threshold = 0.5;
[peaks, peak_locations] = findpeaks(signal, 'MinPeakHeight', threshold);
plot(t, signal);
hold on;
plot(t(peak_locations), peaks, 'ro', 'MarkerSize', 8);
title('Identified Peaks');
xlabel('Time (seconds)');
ylabel('Amplitude');



% Plot the signal with marked cycles

subplot(4, 1, 3);
plot(t, signal);
hold on;
for i = 1:num_cycles
    start_index = cycles(i, 1);
    end_index = cycles(i, 2);
    
    % Ensure indices are within bounds
    start_index = max(1, start_index);
    end_index = min(length(t), end_index);
    
    plot(t(start_index:end_index), signal(start_index:end_index), 'ro-', 'LineWidth', 2);
end
title('Signal with Rainflow Counting Cycles');
xlabel('Time (seconds)');
ylabel('Amplitude');

% Create a 3D histogram
subplot(4, 1, 4);
hist3([t', signal'], [50, 50]); % Adjust the bin numbers as needed
title('3D Histogram of Signal');
xlabel('Time (seconds)');
ylabel('Amplitude');
zlabel('Frequency');

% Adjust subplot spacing
subplot(4, 1, 1);
ylim([-5 5]);
subplot(4, 1, 2);
ylim([-5 5]);
subplot(4, 1, 3);
ylim([-5 5]);

hold off;

%%
% Create a separate figure for the 3D histogram
figure;
hist3([t', signal'], [50, 50]); % Adjust the bin numbers as needed
title('3D Histogram of Signal');
xlabel('Time (seconds)');
ylabel('Amplitude');
zlabel('Frequency');

hold off;
%%
% Create a separate figure for the Rainflow cycles as a 3D scatter plot
figure;
for i = 1:num_cycles
    start_index = cycles(i, 1);
    end_index = cycles(i, 2);
    
    % Ensure indices are within bounds
    start_index = max(1, start_index);
    end_index = min(length(t), end_index);
    
    scatter3(t(start_index), signal(start_index), i, 'filled', 'r');
    hold on;
    scatter3(t(end_index), signal(end_index), i, 'filled', 'b');
end

title('3D Scatter Plot of Rainflow Cycles');
xlabel('Time (seconds)');
ylabel('Amplitude');
zlabel('Cycle Count');
legend('Start of Cycle', 'End of Cycle');

hold off;
