MATLAB Examples

% MATLAB Hello World Examples

% 1. Basic Hello World
disp('Hello, World!');

% 2. Hello World with variable
message = 'Hello, World!';
fprintf('%s\n', message);

% 3. Hello World with function
function result = sayHello()
    result = 'Hello, World!';
end

disp(sayHello());

% 4. Hello World with function parameters
function result = greet(name)
    if nargin < 1
        name = 'World';
    end
    result = sprintf('Hello, %s!', name);
end

disp(greet('World'));
disp(greet('MATLAB'));

% 5. Hello World with formatted output
name = 'MATLAB';
version = 'R2023a';
fprintf('Hello, %s version %s!\n', name, version);

% 6. Hello World multiple times
for i = 1:5
    fprintf('Hello, World! %d\n', i);
end

% 7. Hello World with cell array
greetings = {'Hello', 'Bonjour', 'Hola', 'Ciao', 'こんにちは'};
for i = 1:length(greetings)
    fprintf('%s, World!\n', greetings{i});
end

% 8. Hello World with struct array
translate_greetings.en = 'Hello';
translate_greetings.es = 'Hola';
translate_greetings.fr = 'Bonjour';
translate_greetings.de = 'Hallo';
translate_greetings.ja = 'こんにちは';

fields = fieldnames(translate_greetings);
for i = 1:length(fields)
    lang = fields{i};
    greeting = translate_greetings.(lang);
    fprintf('%s, World! (%s)\n', greeting, lang);
end

% 9. Hello World with condition
hour = hour(datetime('now'));
if hour < 12
    fprintf('Good morning, World!\n');
elseif hour < 18
    fprintf('Good afternoon, World!\n');
else
    fprintf('Good evening, World!\n');
end

% 10. Hello World with user input
% In interactive mode:
% name = input('Enter your name: ');
% fprintf('Hello, %s!\n', name);

% Basic data types examples

% Numbers
integer = 42;
float_num = 3.14159;
double_num = 2.71828;
complex_num = 3 + 4i;

% Strings
text = 'Hello, MATLAB!';
text_with_interpolation = sprintf('The answer is %d', integer);

% Logical
bool_true = true;
bool_false = false;

% Arrays and Matrices
vector = [1, 2, 3, 4, 5];
matrix = [1 2 3; 4 5 6; 7 8 9];

% Cell arrays
cell_array = {1, 'hello', [1 2 3], [1; 2; 3]};

% Struct arrays
struct_array = struct('name', {'Alice', 'Bob'}, 'age', {25, 30});

% Data type demonstrations
fprintf('\n=== Data Types ===\n');
fprintf('Integer: %d, Type: %s\n', integer, class(integer));
fprintf('Float: %f, Type: %s\n', float_num, class(float_num));
fprintf('Complex: %f+%fi, Type: %s\n', real(complex_num), imag(complex_num), class(complex_num));
fprintf('String: %s, Type: %s\n', text, class(text));
fprintf('Boolean: %d, Type: %s\n', logical(bool_true), class(logical(bool_true)));
fprintf('Vector: %s\n', mat2str(vector));
fprintf('Matrix: %s\n', mat2str(matrix));

% Matrix operations
fprintf('\n=== Matrix Operations ===\n');
A = [1 2; 3 4];
B = [5 6; 7 8];

fprintf('Matrix A:\n');
disp(A);
fprintf('Matrix B:\n');
disp(B);

C = A + B;
fprintf('A + B:\n');
disp(C);

D = A * B;
fprintf('A * B:\n');
disp(D);

% Element-wise operations
E = A .* B;
fprintf('A .* B (element-wise):\n');
disp(E);

F = A .^ 2;
fprintf('A .^ 2 (element-wise square):\n');
disp(F);

% Control flow examples
fprintf('\n=== Control Flow ===\n');

% If-else-if-else
age = 18;
if age < 18
    fprintf('You are a minor\n');
elseif age < 65
    fprintf('You are an adult\n');
else
    fprintf('You are a senior\n');
end

% Switch case
day = 3;
switch day
    case 1
        fprintf('Monday\n');
    case 2
        fprintf('Tuesday\n');
    case 3
        fprintf('Wednesday\n');
    case 4
        fprintf('Thursday\n');
    case 5
        fprintf('Friday\n');
    case 6
        fprintf('Saturday\n');
    case 7
        fprintf('Sunday\n');
    otherwise
        fprintf('Invalid day\n');
end

% Loop examples
fprintf('\n=== Loops ===\n');
fprintf('For loop:\n');
fruits = {'apple', 'banana', 'cherry'};
for i = 1:length(fruits)
    fprintf('I like %s\n', fruits{i});
end

fprintf('\nWhile loop:\n');
count = 1;
while count <= 3
    fprintf('Count: %d\n', count);
    count = count + 1;
end

fprintf('\nDo-while loop simulation (not directly supported in MATLAB):\n');
do_count = 1;
while true
    fprintf('Do count: %d\n', do_count);
    if do_count >= 3
        break;
    end
    do_count = do_count + 1;
end

% Function examples
fprintf('\n=== Function Examples ===\n');

% Anonymous functions
add = @(x, y) x + y;
fprintf('5 + 3 = %d\n', add(5, 3));

multiply = @(x, y) x * y;
fprintf('5 * 3 = %d\n', multiply(5, 3));

% Function composition
compose = @(f, g) @(x) f(g(x));
square = @(x) x^2;
double = @(x) 2 * x;
double_square = compose(double, square);
fprintf('Double square of 4 = %d\n', double_square(4));

% String operations
fprintf('\n=== String Operations ===\n');
text = 'Hello, MATLAB Programming!';

fprintf('Original: %s\n', text);
fprintf('Upper case: %s\n', upper(text));
fprintf('Lower case: %s\n', lower(text));
fprintf('Length: %d\n', length(text));
fprintf('Contains MATLAB: %d\n', contains(text, 'MATLAB'));
fprintf('Find MATLAB: %d\n', strfind(text, 'MATLAB'));

% Array operations
fprintf('\n=== Array Operations ===\n');
numbers = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10];

squared = numbers.^2;
fprintf('Squared: %s\n', mat2str(squared));

even_numbers = numbers(mod(numbers, 2) == 0);
fprintf('Even numbers: %s\n', mat2str(even_numbers));

sum_numbers = sum(numbers);
fprintf('Sum of 1-10: %d\n', sum_numbers);

mean_numbers = mean(numbers);
fprintf('Mean of 1-10: %.2f\n', mean_numbers);

% Plotting basics
fprintf('\n=== Basic Plotting ===\n');

% Line plot
figure(1);
x = 0:0.1:2*pi;
y = sin(x);
plot(x, y, 'b-', 'LineWidth', 2);
title('Sine Wave');
xlabel('x');
ylabel('sin(x)');
grid on;

% Multiple lines
figure(2);
x = 0:0.1:2*pi;
y1 = sin(x);
y2 = cos(x);
plot(x, y1, 'b-', x, y2, 'r--', 'LineWidth', 2);
title('Sine and Cosine Waves');
xlabel('x');
ylabel('Value');
legend('sin(x)', 'cos(x)');
grid on;

% Subplots
figure(3);
subplot(2, 2, 1);
plot(x, sin(x));
title('Sin(x)');

subplot(2, 2, 2);
plot(x, cos(x));
title('Cos(x)');

subplot(2, 2, 3);
plot(x, tan(x));
title('Tan(x)');

subplot(2, 2, 4);
plot(x, exp(-x));
title('e^(-x)');
sgtitle('Trigonometric and Exponential Functions');

% 3D Plotting
fprintf('\n=== 3D Plotting ===\n');
figure(4);
[X, Y] = meshgrid(-2:0.2:2, -2:0.2:2);
Z = X .* exp(-X.^2 - Y.^2);
surf(X, Y, Z);
title('3D Surface Plot');
xlabel('X');
ylabel('Y');
zlabel('Z');
colorbar;

% Working with Files
fprintf('\n=== File Operations ===\n');

% Write to file
filename = 'matlab_hello.txt';
fid = fopen(filename, 'w');
if fid > 0
    fprintf(fid, 'Hello from MATLAB!\n');
    fprintf(fid, 'This is line 2\n');
    fprintf(fid, 'Current time: %s\n', datestr(now));
    fclose(fid);
    fprintf('File written successfully: %s\n', filename);
else
    fprintf('Failed to create file: %s\n', filename);
end

% Read from file
if exist(filename, 'file')
    fid = fopen(filename, 'r');
    if fid > 0
        fprintf('Reading file contents:\n');
        while ~feof(fid)
            line = fgetl(fid);
            if ~ischar(line)  % Handle empty lines
                break;
            end
            fprintf('%s\n', line);
        end
        fclose(fid);
    end
else
    fprintf('File does not exist: %s\n', filename);
end

% Error handling
fprintf('\n=== Error Handling ===\n');

% Try-catch pattern (using try-catch)
try
    result = 10 / 0;
    fprintf('Division result: %f\n', result);
catch ME
    fprintf('Error caught: %s\n', ME.message);
end

% Input validation
function validated_divide(a, b)
    if b == 0
        error('Division by zero is not allowed!');
    end
    result = a / b;
    fprintf('%g / %g = %g\n', a, b, result);
end

try
    validated_divide(10, 2);
    validated_divide(10, 0);
catch ME
    fprintf('Validation error: %s\n', ME.message);
end

% Data structures
fprintf('\n=== Data Structures ===\n');

% Cell array operations
cell_array = {1, 'hello', [1 2 3], struct('field1', 'value1')};
fprintf('Cell array size: %s\n', mat2str(size(cell_array)));
for i = 1:length(cell_array)
    fprintf('Cell %d: %s\n', i, class(cell_array{i}));
end

% Struct array operations
people = struct('name', {}, 'age', {});
people(1) = struct('name', 'Alice', 'age', 25);
people(2) = struct('name', 'Bob', 'age', 30);
people(3) = struct('name', 'Charlie', 'age', 35);

fprintf('People array:\n');
for i = 1:length(people)
    fprintf('%s: %d years old\n', people(i).name, people(i).age);
end

% Time and date operations
fprintf('\n=== Time and Date Operations ===\n');
current_time = now;
fprintf('Current time: %s\n', datestr(current_time));
fprintf('Year: %d\n', year(current_time));
fprintf('Month: %d\n', month(current_time));
fprintf('Day: %d\n', day(current_time));
fprintf('Hour: %d\n', hour(current_time));
fprintf('Minute: %d\n', minute(current_time));
fprintf('Second: %d\n', second(current_time));

% Time calculations
start_time = tic;
pause(2);  % Wait 2 seconds
elapsed_time = toc(start_time);
fprintf('Elapsed time: %.2f seconds\n', elapsed_time);

% Built-in functions
fprintf('\n=== Built-in Functions ===\n');

% Mathematical functions
numbers = [1, 2, 3, 4, 5];
fprintf('Numbers: %s\n', mat2str(numbers));
fprintf('Sum: %d\n', sum(numbers));
fprintf('Mean: %.2f\n', mean(numbers));
fprintf('Standard deviation: %.2f\n', std(numbers));
fprintf('Maximum: %d\n', max(numbers));
fprintf('Minimum: %d\n', min(numbers));

% Statistical functions
data = [10, 15, 12, 18, 20, 14, 16, 13, 19, 11];
fprintf('\nStatistical analysis of data:\n');
fprintf('Mean: %.2f\n', mean(data));
fprintf('Median: %.2f\n', median(data));
fprintf('Mode: %s\n', mat2str(mode(data)));
fprintf('Range: %d\n', max(data) - min(data));
fprintf('Variance: %.2f\n', var(data));
fprintf('Standard deviation: %.2f\n', std(data));

fprintf('\n=== All MATLAB Hello World Examples Complete! ===');

% MATLAB Numerical Analysis Examples

% 1. Numerical Integration and Differentiation

% Numerical integration using trapezoidal rule
function I = trapezoidal_rule(f, a, b, n)
    % f: function handle
    % a, b: integration limits
    % n: number of subintervals

    h = (b - a) / n;
    x = a:h:b;
    y = arrayfun(f, x);

    I = h * (0.5 * y(1) + sum(y(2:end-1)) + 0.5 * y(end));
end

% Numerical differentiation using forward difference
function df = forward_difference(f, x, h)
    df = (f(x + h) - f(x)) / h;
end

% Central difference (more accurate)
function df = central_difference(f, x, h)
    df = (f(x + h) - f(x - h)) / (2 * h);
end

% Test integration and differentiation
f = @(x) x.^2;
a = 0;
b = 2;
n = 100;

integral_result = trapezoidal_rule(f, a, b, n);
true_integral = (b^3 - a^3) / 3;

fprintf('Numerical integration of x^2 from 0 to 2:\n');
fprintf('Numerical result: %.6f\n', integral_result);
fprintf('True result: %.6f\n', true_integral);
fprintf('Error: %.6f\n', abs(integral_result - true_integral));

x0 = 1;
h = 0.001;
f_prime_exact = 2 * x0;
f_prime_forward = forward_difference(f, x0, h);
f_prime_central = central_difference(f, x0, h);

fprintf('\nDerivative of x^2 at x=1:\n');
fprintf('Exact derivative: %.6f\n', f_prime_exact);
fprintf('Forward difference: %.6f\n', f_prime_forward);
fprintf('Central difference: %.6f\n', f_prime_central);
fprintf('Forward error: %.6f\n', abs(f_prime_forward - f_prime_exact));
fprintf('Central error: %.6f\n', abs(f_prime_central - f_prime_exact));

% 2. Root Finding Methods

% Bisection method
function root = bisection_method(f, a, b, tol, max_iter)
    if f(a) * f(b) >= 0
        error('Function must have opposite signs at endpoints');
    end

    iter = 0;
    while abs(b - a) > tol && iter < max_iter
        c = (a + b) / 2;
        if f(c) == 0
            break;
        elseif f(a) * f(c) < 0
            b = c;
        else
            a = c;
        end
        iter = iter + 1;
    end

    root = (a + b) / 2;
    fprintf('Bisection method converged in %d iterations\n', iter);
end

% Newton's method
function [root, iterations] = newton_method(f, df, x0, tol, max_iter)
    x = x0;
    iter = 0;

    while iter < max_iter
        fx = f(x);
        dfx = df(x);

        if abs(dfx) < tol
            warning('Derivative too small');
            break;
        end

        x_new = x - fx / dfx;

        if abs(x_new - x) < tol
            break;
        end

        x = x_new;
        iter = iter + 1;
    end

    root = x;
    iterations = iter;
end

% Test root finding
f_root = @(x) x^3 - 2*x - 5;
df_root = @(x) 3*x^2 - 2;

fprintf('\n=== Root Finding Methods ===\n');

% Bisection method
fprintf('Bisection method for f(x) = x^3 - 2x - 5\n');
root_bisection = bisection_method(f_root, 2, 3, 1e-6, 100);
fprintf('Root: %.6f\n', root_bisection);
fprintf('f(root): %.6f\n', f_root(root_bisection));

% Newton's method
fprintf('\nNewton''s method for f(x) = x^3 - 2x - 5\n');
[root_newton, iterations] = newton_method(f_root, df_root, 2, 1e-6, 100);
fprintf('Root: %.6f\n', root_newton);
fprintf('f(root): %.6f\n', f_root(root_newton));
fprintf('Iterations: %d\n', iterations);

% 3. Linear Systems

% Gaussian elimination
function x = gaussian_elimination(A, b)
    [n, m] = size(A);
    [n_b, m_b] = size(b);

    if n ~= m || n ~= n_b
        error('Matrix must be square and compatible with b');
    end

    % Forward elimination
    for k = 1:n-1
        for i = k+1:n
            factor = A(i,k) / A(k,k);
            A(i,k:end) = A(i,k:end) - factor * A(k,k:end);
            b(i) = b(i) - factor * b(k);
        end
    end

    % Back substitution
    x = zeros(n, 1);
    for i = n:-1:1
        x(i) = (b(i) - A(i,i+1:n) * x(i+1:n)) / A(i,i);
    end
end

% Test linear system solving
fprintf('\n=== Linear System Solving ===\n');

A = [3 1 2; 1 4 1; 2 2 3];
b = [7; 6; 8];

fprintf('Coefficient matrix A:\n');
disp(A);
fprintf('Right hand side b:\n');
disp(b);

% Gaussian elimination
x_ge = gaussian_elimination(A, b);
fprintf('Solution by Gaussian elimination:\n');
disp(x_ge);

% MATLAB built-in solver
x_builtin = A \ b;
fprintf('Solution by MATLAB solver:\n');
disp(x_builtin);

% Verification
fprintf('Verification (A*x):\n');
disp(A * x_ge);

% 4. Interpolation

% Lagrange interpolation
function y_interp = lagrange_interpolation(x_data, y_data, x_interp)
    n = length(x_data);
    y_interp = 0;

    for i = 1:n
        L_i = 1;
        for j = 1:n
            if i ~= j
                L_i = L_i * (x_interp - x_data(j)) / (x_data(i) - x_data(j));
            end
        end
        y_interp = y_interp + y_data(i) * L_i;
    end
end

% Test interpolation
fprintf('\n=== Interpolation ===\n');

x_data = [0, 1, 2, 3, 4];
y_data = [0, 1, 8, 27, 64];  % y = x^3
x_interp = 2.5;

y_interp = lagrange_interpolation(x_data, y_data, x_interp);
fprintf('Interpolation at x = %.1f\n', x_interp);
fprintf('Interpolated y = %.6f\n', y_interp);
fprintf('True y = %.6f\n', x_interp^3);

% Plot interpolation
figure(1);
x_plot = linspace(0, 4, 100);
y_interp_plot = zeros(size(x_plot));
for i = 1:length(x_plot)
    y_interp_plot(i) = lagrange_interpolation(x_data, y_data, x_plot(i));
end

plot(x_data, y_data, 'ro', 'MarkerSize', 8, 'LineWidth', 2);
hold on;
plot(x_plot, y_interp_plot, 'b-', 'LineWidth', 1);
plot(x_interp, y_interp, 'ks', 'MarkerSize', 8, 'LineWidth', 2);
title('Lagrange Interpolation');
xlabel('x');
ylabel('y');
legend('Data Points', 'Interpolated Curve', 'Interpolated Point');
grid on;
hold off;

% 5. Curve Fitting

% Linear regression
function [slope, intercept] = linear_regression(x, y)
    n = length(x);
    sum_x = sum(x);
    sum_y = sum(y);
    sum_xy = sum(x .* y);
    sum_x2 = sum(x .^2);

    slope = (n * sum_xy - sum_x * sum_y) / (n * sum_x2 - sum_x^2);
    intercept = (sum_y - slope * sum_x) / n;
end

% Polynomial fitting
function coeffs = polynomial_fit(x, y, degree)
    coeffs = polyfit(x, y, degree);
end

% Test curve fitting
fprintf('\n=== Curve Fitting ===\n');

x_fit = [0; 1; 2; 3; 4; 5];
y_fit = [1.1; 1.9; 3.1; 4.0; 5.1; 5.9];  % Approximately y = x + 1

% Linear regression
[slope, intercept] = linear_regression(x_fit, y_fit);
fprintf('Linear regression: y = %.3fx + %.3f\n', slope, intercept);

% Polynomial fitting (quadratic)
poly_coeffs = polynomial_fit(x_fit, y_fit, 2);
fprintf('Quadratic fit: y = %.3fx^2 + %.3fx + %.3f\n', poly_coeffs(1), poly_coeffs(2), poly_coeffs(3));

% Plot fitting results
figure(2);
plot(x_fit, y_fit, 'o', 'MarkerSize', 8, 'LineWidth', 2);
hold on;

% Linear fit
x_lin = linspace(0, 5, 100);
y_lin = slope * x_lin + intercept;
plot(x_lin, y_lin, 'r-', 'LineWidth', 2);

% Polynomial fit
x_poly = linspace(0, 5, 100);
y_poly = polyval(poly_coeffs, x_poly);
plot(x_poly, y_poly, 'b-', 'LineWidth', 2);

title('Curve Fitting Comparison');
xlabel('x');
ylabel('y');
legend('Data Points', 'Linear Fit', 'Quadratic Fit');
grid on;
hold off;

% 6. Ordinary Differential Equations

% Euler's method
function [t, y] = euler_method(f, t0, y0, h, n)
    t = zeros(n+1, 1);
    y = zeros(n+1, 1);

    t(1) = t0;
    y(1) = y0;

    for i = 1:n
        y(i+1) = y(i) + h * f(t(i), y(i));
        t(i+1) = t(i) + h;
    end
end

% Runge-Kutta 4th order
function [t, y] = rk4(f, t0, y0, h, n)
    t = zeros(n+1, 1);
    y = zeros(n+1, 1);

    t(1) = t0;
    y(1) = y0;

    for i = 1:n
        k1 = f(t(i), y(i));
        k2 = f(t(i) + h/2, y(i) + h*k1/2);
        k3 = f(t(i) + h/2, y(i) + h*k2/2);
        k4 = f(t(i) + h, y(i) + h*k3);

        y(i+1) = y(i) + h * (k1 + 2*k2 + 2*k3 + k4) / 6;
        t(i+1) = t(i) + h;
    end
end

% Test ODE solvers
fprintf('\n=== ODE Solvers ===\n');

% dy/dt = -2*y + 1, y(0) = 0
f_ode = @(t, y) -2*y + 1;
t0 = 0;
y0 = 0;
h = 0.1;
n = 50;

% Euler method
[t_euler, y_euler] = euler_method(f_ode, t0, y0, h, n);

% Runge-Kutta method
[t_rk4, y_rk4] = rk4(f_ode, t0, y0, h, n);

% Exact solution
y_exact = @(t) 0.5 - 0.5 * exp(-2*t);

% Plot results
figure(3);
plot(t_euler, y_euler, 'r-', 'LineWidth', 2);
hold on;
plot(t_rk4, y_rk4, 'b-', 'LineWidth', 2);
plot(t_euler, y_exact(t_euler), 'k--', 'LineWidth', 2);
title('ODE Solution Comparison');
xlabel('t');
ylabel('y');
legend('Euler Method', 'RK4 Method', 'Exact Solution');
grid on;
hold off;

% Error analysis
error_euler = abs(y_euler - y_exact(t_euler));
error_rk4 = abs(y_rk4 - y_exact(t_rk4));

fprintf('Maximum error (Euler): %.6f\n', max(error_euler));
fprintf('Maximum error (RK4): %.6f\n', max(error_rk4));

% 7. Fourier Analysis

% Discrete Fourier Transform
function X = dft_custom(x)
    N = length(x);
    X = zeros(N, 1);

    for k = 0:N-1
        sum_val = 0;
        for n = 0:N-1
            sum_val = sum_val + x(n+1) * exp(-1i * 2 * pi * k * n / N);
        end
        X(k+1) = sum_val;
    end
end

% Signal processing example
fprintf('\n=== Fourier Analysis ===\n');

% Create test signal
fs = 1000;  % Sampling frequency
t = 0:1/fs:1-1/fs*100;  % 1 second of data
f1 = 50;    % 50 Hz
f2 = 120;   % 120 Hz
signal = sin(2*pi*f1*t) + 0.5*sin(2*pi*f2*t) + 0.2*randn(size(t));

% Compute FFT
X = fft(signal);
n = length(signal);
freq = (0:n-1)*(fs/n);

% Plot signal and spectrum
figure(4);
subplot(2, 1, 1);
plot(t, signal);
title('Original Signal');
xlabel('Time (s)');
ylabel('Amplitude');
grid on;

subplot(2, 1, 2);
plot(freq(1:n/2), abs(X(1:n/2)));
title('Single-Sided Amplitude Spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
grid on;

% 8. Optimization

% Golden section search for 1D optimization
function [x_min, f_min] = golden_section(f, a, b, tol)
    phi = (sqrt(5) - 1) / 2;

    while abs(b - a) > tol
        c = b - phi * (b - a);
        d = a + phi * (b - a);

        if f(c) < f(d)
            b = d;
        else
            a = c;
        end
    end

    x_min = (a + b) / 2;
    f_min = f(x_min);
end

% Test optimization
fprintf('\n=== Optimization ===\n');

% Minimize f(x) = x^2 - 4x + 3 = (x-1)(x-3)
f_opt = @(x) x^2 - 4*x + 3;
a_opt = 0;
b_opt = 5;

[x_min, f_min] = golden_section(f_opt, a_opt, b_opt, 1e-6);
fprintf('Minimum found at x = %.6f, f(x) = %.6f\n', x_min, f_min);
fprintf('True minimum at x = %.1f, f(x) = %.1f\n', 1.0, 1.0);

% Plot optimization function
figure(5);
x_plot = a_opt:0.01:b_opt;
y_plot = arrayfun(f_opt, x_plot);
plot(x_plot, y_plot, 'b-', 'LineWidth', 2);
hold on;
plot(x_min, f_min, 'ro', 'MarkerSize', 10, 'LineWidth', 2);
title('One-D Optimization');
xlabel('x');
ylabel('f(x)');
legend('Function', 'Minimum');
grid on;
hold off;

% 9. Random Number Generation and Monte Carlo

% Monte Carlo estimation of pi
function pi_estimate = monte_carlo_pi(n)
    points_inside = 0;

    for i = 1:n
        x = rand();  % Uniform [0,1)
        y = rand();  % Uniform [0,1]
        if x^2 + y^2 <= 1
            points_inside = points_inside + 1;
        end
    end

    pi_estimate = 4 * points_inside / n;
end

% Test Monte Carlo
fprintf('\n=== Monte Carlo Methods ===\n');

n_trials = [100, 1000, 10000, 100000];
for n = n_trials
    pi_est = monte_carlo_pi(n);
    error = abs(pi_est - pi);
    fprintf('N = %d: π ≈ %.6f (error = %.6f)\n', n, pi_est, error);
end

% 10. Numerical Stability

% Condition number estimation
function cond_num = estimate_condition_number(A)
    try
        [V, D] = eig(A);
        singular_values = abs(diag(D));
        if min(singular_values) > 0
            cond_num = max(singular_values) / min(singular_values);
        else
            cond_num = Inf;
        end
    catch
        cond_num = Inf;
    end
end

% Test condition number
fprintf('\n=== Numerical Stability ===\n');

% Well-conditioned matrix
A_well = [1 2; 3 4];
cond_well = cond(A_well);

% Ill-conditioned matrix
A_ill = [1 2; 1.0000001 2];
cond_ill = cond(A_ill);

fprintf('Well-conditioned matrix condition number: %.2f\n', cond_well);
fprintf('Ill-conditioned matrix condition number: %.2e\n', cond_ill);

fprintf('\n=== All Numerical Analysis Examples Complete! ===');

% MATLAB Signal Processing Examples

% 1. Digital Filter Design and Analysis

% FIR Filter Design using Window Method
function h = fir_filter_design(N, cutoff, fs, window_type)
    % N: filter order
    % cutoff: cutoff frequency (Hz)
    % fs: sampling frequency (Hz)
    % window_type: 'rectangular', 'hamming', 'hanning', 'blackman'

    normalized_cutoff = cutoff / (fs/2);

    % Ideal impulse response
    h_ideal = zeros(N, 1);
    for n = 0:N-1
        if n ~= N/2
            h_ideal(n+1) = 2 * normalized_cutoff * ...
                sinc(2 * normalized_cutoff * (n - N/2));
        end
    end

    % Apply window
    switch lower(window_type)
        case 'rectangular'
            w = ones(N, 1);
        case 'hamming'
            w = hamming(N);
        case 'hanning'
            w = hanning(N);
        case 'blackman'
            w = blackman(N);
        otherwise
            error('Unknown window type');
    end

    h = h_ideal .* w';
    h = h / sum(h);  % Normalize
end

% Filter frequency response
function [H, freq] = freq_response(h, fs, nfft)
    H = fft(h, nfft);
    freq = (0:nfft-1) * (fs/nfft);
end

% Test FIR filter design
fprintf('\n=== FIR Filter Design ===\n');

fs = 1000;  % Sampling frequency
cutoff = 200;  % Cutoff frequency
N = 51;      % Filter order

h_hamming = fir_filter_design(N, cutoff, fs, 'hamming');
[H, freq] = freq_response(h_hamming, fs, 1024);

% Plot frequency response
figure(1);
subplot(2, 1, 1);
plot(freq(1:512), abs(H(1:512)), 'b-', 'LineWidth', 2);
title('Hamming Window FIR Filter - Magnitude Response');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
grid on;

subplot(2, 1, 2);
plot(freq(1:512), 20*log10(abs(H(1:512))), 'r-', 'LineWidth', 2);
title('Hamming Window FIR Filter - Magnitude Response (dB)');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
grid on;

% Impulse response
figure(2);
stem(0:N-1, h_hamming, 'filled');
title('Hamming Window FIR Filter - Impulse Response');
xlabel('Sample');
ylabel('Amplitude');
grid on;

% 2. IIR Filter Design

% Butterworth filter design
function [b, a] = butterworth_filter(order, cutoff, fs, filter_type)
    normalized_cutoff = cutoff / (fs/2);

    if strcmp(filter_type, 'low') == 1
        [b, a] = butter(order, normalized_cutoff, 'low');
    elseif strcmp(filter_type, 'high') == 1
        [b, a] = butter(order, normalized_cutoff, 'high');
    elseif strcmp(filter_type, 'band') == 1
        % Bandpass would require cutoffs array
        error('Bandpass not implemented in this simple version');
    else
        error('Unknown filter type');
    end
end

% Test IIR filter
fprintf('\n=== IIR Filter Design ===\n');

order = 4;
cutoff_iir = 150;
[b, a] = butterworth_filter(order, cutoff_iir, fs, 'low');

% Frequency response
[H_iir, freq_iir] = freqz(b, a, 1024, fs);

figure(3);
subplot(2, 1, 1);
plot(freq_iir(1:512), abs(H_iir(1:512)), 'g-', 'LineWidth', 2);
title('Butterworth IIR Filter - Magnitude Response');
xlabel('Frequency (Hz)');
ylabel('Magnitude');
grid on;

subplot(2, 1, 2);
plot(freq_iir(1:512), 20*log10(abs(H_iir(1:512))), 'r-', 'LineWidth', 2);
title('Butterworth IIR Filter - Magnitude Response (dB)');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
grid on;

% Pole-zero plot
figure(4);
zplane(b, a);
title('Pole-Zero Plot');
legend('Poles', 'Zeros');
grid on;

% 3. Convolution and Correlation

% Convolution example
fprintf('\n=== Convolution and Correlation ===\n');

% Create test signals
x = [1, 2, 3, 4, 5];
h = [1, 1, 1];  % Simple moving average filter

% Convolution
y = conv(x, h);
fprintf('Input signal: %s\n', mat2str(x));
fprintf('Filter impulse response: %s\n', mat2str(h));
fprintf('Convolution result: %s\n', mat2str(y));

% Circular convolution
circular_y = cconv(x, h, length(x) + length(h) - 1);
fprintf('Circular convolution result: %s\n', mat2str(circular_y));

% Cross-correlation
cross_corr = xcorr(x, x);
fprintf('Autocorrelation of x: %s\n', mat2str(cross_corr(length(cross_corr)/2+1:end)));

% Plot convolution results
figure(5);
subplot(3, 1, 1);
stem(x, 'b-o', 'LineWidth', 2);
title('Input Signal');
xlabel('Sample');
ylabel('Amplitude');
grid on;

subplot(3, 1, 2);
stem(h, 'r-o', 'LineWidth', 2);
title('Filter Impulse Response');
xlabel('Sample');
ylabel('Amplitude');
grid on;

subplot(3, 1, 3);
stem(y, 'g-o', 'LineWidth', 2);
title('Convolution Result');
xlabel('Sample');
ylabel('Amplitude');
grid on;

% 4. Spectral Analysis

% Power Spectral Density estimation using Welch's method
function psd = welch_psd(x, fs, nfft, noverlap)
    n = length(x);
    if n < nfft
        x = [x, zeros(1, nfft-n)];
        n = nfft;
    end

    % Divide into segments
    segment_length = floor(nfft);
    nsegments = floor((length(x) - noverlap) / (segment_length - noverlap));

    psd = zeros(nfft, 1);

    for i = 1:nsegments
        start_idx = (i-1) * (segment_length - noverlap) + 1;
        end_idx = start_idx + segment_length - 1;
        segment = x(start_idx:end_idx);

        % Apply window (Hamming)
        window = hamming(length(segment));
        segment = segment .* window;

        % FFT and power calculation
        X = fft(segment, nfft);
        power = abs(X).^2 / (fs * sum(window.^2));

        psd = psd + power;
    end

    psd = psd / nsegments;
    psd = [psd(1:nfft/2); psd(nfft/2:-1:1)];  % Make single-sided
end

% Generate test signal with multiple frequencies
fprintf('\n=== Spectral Analysis ===\n');

fs_psd = 2000;
t_psd = 0:1/fs_psd:1-1/fs_psd*2;  % 2 seconds
signal_psd = sin(2*pi*100*t_psd) + 0.5*sin(2*pi*250*t_psd) + ...
              0.3*sin(2*pi*400*t_psd) + 0.1*randn(size(t_psd));

% Compute PSD
psd = welch_psd(signal_psd, fs_psd, 1024, 512);
freq_psd = (0:length(psd)-1) * (fs_psd/2) / length(psd);

figure(6);
plot(freq_psd/1000, psd, 'b-', 'LineWidth', 2);
title('Power Spectral Density (Welch''s Method)');
xlabel('Frequency (kHz)');
ylabel('Power/Frequency (dB/Hz)');
grid on;

% 5. Window Functions and Spectral Leakage

fprintf('\n=== Window Functions Analysis ===\n');

n_fft = 256;
freq_fft = (0:n_fft-1) * (fs/2) / n_fft;

% Test different windows
test_signal = sin(2*pi*100*t_psd(1:n_fft));

windows = {'rectangular', 'hamming', 'hanning', 'blackman'};

figure(7);
for i = 1:length(windows)
    window_name = windows{i};

    switch window_name
        case 'rectangular'
            w = rectwin(n_fft);
        case 'hamming'
            w = hamming(n_fft);
        case 'hanning'
            w = hann(n_fft);
        case 'blackman'
            w = blackman(n_fft);
    end

    windowed_signal = test_signal .* w;
    X = fft(windowed_signal);

    subplot(2, 2, i);
    plot(freq_fft/1000, 20*log10(abs(X(1:n_fft/2))), 'LineWidth', 2);
    title([window_name ' Window - Spectrum']);
    xlabel('Frequency (kHz)');
    ylabel('Magnitude (dB)');
    grid on;
end

sgtitle('Window Function Comparison');

% 6. Digital Communications Simulation

% ASK modulation
function [ask_signal, carrier] = ask_modulation(data_bits, fs, bit_rate, carrier_freq)
    % data_bits: binary data
    fs: sampling frequency
    bit_rate: bits per second
    carrier_freq: carrier frequency

    samples_per_bit = fs / bit_rate;
    t_total = length(data_bits) * samples_per_bit / fs;
    t = 0:1/fs:t_total-1;

    % Generate NRZ signal from bits
    nrz_signal = zeros(size(t));
    for i = 1:length(data_bits)
        start_idx = (i-1) * samples_per_bit + 1;
        end_idx = i * samples_per_bit;
        nrz_signal(start_idx:end_idx) = data_bits(i);
    end

    % ASK modulation
    carrier = cos(2*pi*carrier_freq*t);
    ask_signal = nrz_signal .* (2*carrier - 1);
end

% Generate test data and modulate
fprintf('\n=== Digital Communications ===\n');

data_bits = [1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0];
fs_comm = 10000;  % 10 kHz
bit_rate = 1000;  % 1 kbps
carrier_freq = 2000; % 2 kHz

[ask_signal, carrier] = ask_modulation(data_bits, fs_comm, bit_rate, carrier_freq);

% Plot ASK signal
figure(8);
subplot(2, 1, 1);
plot(0:length(ask_signal)/fs_comm-1/fs_comm, ask_signal, 'b-', 'LineWidth', 1);
title('ASK Modulated Signal');
xlabel('Time (s)');
ylabel('Amplitude');
grid on;

subplot(2, 1, 2);
% Plot spectrum of ASK signal
X_ask = fft(ask_signal);
freq_ask = (0:length(X_ask)-1) * (fs_comm/2) / length(X_ask);
plot(freq_ask(1:length(X_ask)/2)/1000, abs(X_ask(1:length(X_ask)/2)), 'r-', 'LineWidth', 1);
title('ASK Signal Spectrum');
xlabel('Frequency (kHz)');
ylabel('Magnitude');
grid on;

% 7. Image Processing

% Simple image operations
fprintf('\n=== Image Processing ===\n');

% Create a test image (sine wave pattern)
image_size = 256;
[X, Y] = meshgrid(1:image_size, 1:image_size);
test_image = sin(2*pi*X/50) .* cos(2*pi*Y/50);

% Add noise
noisy_image = test_image + 0.1 * randn(image_size);

% Apply filter (simple averaging filter)
filter_kernel = ones(5, 5) / 25;
filtered_image = imfilter(noisy_image, filter_kernel);

% Calculate SNR
signal_power = mean(test_image(:).^2);
noise_power = mean((noisy_image - test_image).^2);
snr_db = 10 * log10(signal_power / noise_power);

fprintf('Original SNR: %.2f dB\n', snr_db);

signal_power_filtered = mean(test_image(:).^2);
noise_power_filtered = mean((filtered_image - test_image).^2);
snr_filtered_db = 10 * log10(signal_power_filtered / noise_power_filtered);

fprintf('Filtered SNR: %.2f dB\n', snr_filtered_db);

% Display images
figure(9);
subplot(2, 2, 1);
imshow(test_image, []);
title('Original Image');
axis off;

subplot(2, 2, 2);
imshow(noisy_image, []);
title('Noisy Image');
axis off;

subplot(2, 2, 3);
imshow(filtered_image, []);
title('Filtered Image');
axis off;

% Display difference
subplot(2, 2, 4);
imshow(abs(test_image - filtered_image), []);
title('Difference Image');
axis off;

sgtitle('Image Filtering Example');

% 8. Adaptive Filtering

% LMS adaptive filter for noise cancellation
function [e, w] = lms_adaptive_filter(x, d, mu, n_taps)
    % x: reference signal
    % d: desired signal
    % mu: step size
    % n_taps: number of filter taps

    w = zeros(n_taps, 1);  % Initial weights
    e = zeros(length(x), 1);  % Error signal

    % Buffer for input
    x_buffer = zeros(n_taps, 1);

    for i = n_taps:length(x)
        % Update buffer
        x_buffer = [x(i); x_buffer(1:end-1)];

        % Filter output
        y = w' * x_buffer;

        % Error signal
        e(i) = d(i) - y;

        % Update weights (LMS algorithm)
        w = w + mu * e(i) * x_buffer;
    end
end

% Test adaptive filtering
fprintf('\n=== Adaptive Filtering ===\n');

% Create reference noise and desired signal
rng(42); % Set random seed for reproducibility
noise = 0.5 * randn(1000, 1);
clean_signal = sin(2*pi*50*(1:1000)/1000);
desired_signal = clean_signal(1:1000)';
reference_noise = [0; noise(1:999)]; % Pre-pend zero for alignment

% LMS parameters
mu = 0.01;      % Step size
n_taps = 20;     % Number of filter taps

% Apply adaptive filter
[e, w] = lms_adaptive_filter(reference_noise, desired_signal, mu, n_taps);

% Plot results
figure(10);
subplot(3, 1, 1);
plot(desired_signal, 'b-', 'LineWidth', 1);
title('Desired Signal');
xlabel('Sample');
ylabel('Amplitude');
grid on;

subplot(3, 1, 2);
plot(clean_signal + noise, 'r-', 'LineWidth', 1);
hold on;
plot(e, 'g-', 'LineWidth', 1);
title('Noise and Error Signal');
xlabel('Sample');
ylabel('Amplitude');
legend('Noisy Signal', 'Error', 'Location', 'northwest');
grid on;
hold off;

subplot(3, 1, 3);
plot(w, 'm-', 'LineWidth', 2);
title('Adaptive Filter Weights');
xlabel('Tap Index');
ylabel('Weight Value');
grid on;

% Calculate final SNR
final_snr = 10 * log10(var(desired_signal) / var(e));
fprintf('Final SNR after filtering: %.2f dB\n', final_snr);

% 9. Filter Banks and Multirate Systems

% Create an analysis filter bank
fprintf('\n=== Filter Bank ===\n');

% Design analysis filters (octave band filters)
center_freqs = [100, 200, 400, 800, 1600];
bandwidths = center_freqs * [0.7, 0.7, 0.7, 0.7, 0.7];

n_bands = length(center_freqs);
filters = cell(n_bands, 1);

for i = 1:n_bands
    % Design bandpass filter for each band
    [b, a] = butterworth_filter(4, ...
        [center_freqs(i) - bandwidths(i)/2, center_freqs(i) + bandwidths(i)/2], ...
        fs_comm, 'band');
    filters{i} = struct('b', b, 'a', a, 'center_freq', center_freq(i), 'bandwidth', bandwidths(i));
end

% Test signal with multiple frequency components
multitone_signal = sin(2*pi*100*t_psd) + ...
                    0.5*sin(2*pi*200*t_psd) + ...
                    0.3*sin(2*pi*400*t_psd) + ...
                    0.2*sin(2*pi*800*t_psd);

% Apply filter bank
subband_signals = cell(n_bands, 1);
figure(11);

for i = 1:n_bands
    subband_signals{i} = filter(filters{i}.b, filters{i}.a, multitone_signal);
end

% Plot subband signals
for i = 1:n_bands
    subplot(3, 2, i);
    plot(0:length(subband_signals{i})/fs_comm-1/fs_comm, subband_signals{i}, 'LineWidth', 1);
    title(sprintf('Subband %d (%.0f Hz)', i, center_freqs(i)));
    xlabel('Time (s)');
    ylabel('Amplitude');
    grid on;
end

sgtitle('Filter Bank Output');

% 10. Real-Time Signal Processing Simulation

% Simple real-time processing simulation
fprintf('\n=== Real-Time Signal Processing Simulation ===\n');

buffer_size = 256;
overlap = 128;
hop_size = buffer_size - overlap;

% Initialize buffer and windows
buffer = zeros(buffer_size, 1);
window = hamming(buffer_size);
n_windows = 0;

% Process signal in overlapping windows
total_samples = length(signal_psd);
processed_frames = zeros(buffer_size, ceil((total_samples - buffer_size) / hop_size) + 1);

for frame_idx = 1:ceil((total_samples - buffer_size) / hop_size) + 1
    start_idx = (frame_idx - 1) * hop_size + 1;
    end_idx = min(start_idx + buffer_size - 1, total_samples);

    if end_idx > buffer_size
        % Overlap-add method
        buffer(1:buffer_size-start_idx) = buffer(start_idx:end_idx);

        % Process window
        windowed_buffer = buffer .* window;
        processed_frame = fft(windowed_buffer);

        % Simple processing: frequency thresholding
        processed_frame(abs(processed_frame) > 5) = 0;

        if n_windows > 0
            % Overlap-add to previous frame
            processed_frames(:, n_windows) = processed_frames(:, n_windows) + processed_frame(1:buffer_size);
        else
            processed_frames(:, 1) = processed_frame(1:buffer_size);
        end

        n_windows = n_windows + 1;

        buffer = zeros(buffer_size, 1);
    else
        buffer(1:end_idx) = signal_psd(start_idx:end_idx);
    end
end

fprintf('Processed %d frames\n', n_windows);

% Plot processed frames
figure(12);
if n_windows > 0
    first_frame = processed_frames(:, 1);
    last_frame = processed_frames(:, n_windows);

    subplot(2, 1, 1);
    plot(0:length(first_frame)/fs_comm-1/fs_comm, first_frame, 'b-', 'LineWidth', 1);
    title('First Processed Frame');
    xlabel('Time (s)');
    ylabel('Amplitude');
    grid on;

    subplot(2, 1, 2);
    plot(0:length(last_frame)/fs_comm-1/fs_comm, last_frame, 'r-', 'LineWidth', 1);
    title('Last Processed Frame');
    xlabel('Time (s)');
    ylabel('Amplitude');
    grid on;
end

sgtitle('Real-Time Processing Results');

fprintf('\n=== All Signal Processing Examples Complete! ===');