MATLAB: 循环loop、脚本-文件和绘图、函数-文件（M-文件）、信号生成和卷积、傅里叶变换、滤波器filter的实现

引言

MATLAB软件具备多种功能，在数据分析方面，它可以对数据进行探查、建模和可视化；在创建图形方面，MATLAB可以通过图形可视化并探查数据；在编程方面，MATLAB可以进行编程创建脚本、函数和类。同时，MATLAB可以与多种编程语言使用，例如可以将 MATLAB 与 Python、C/C++、Fortran、Java 及其他语言结合使用；在App 构建方面，MATLAB可以创建桌面 App 和 Web App；在硬件方面，MATLAB可以连接到硬件；在并行计算方面，MATLAB使用多核台式机、GPU、集群、网格和云来进行大规模计算...

 

1.0 Loops循环

1）通过单击文件→新建→空白M-文件打开新脚本。（Open a new script by clicking File → New → Blank M-File.）

2）设置MATLAB的路径以运行文件（文件→设置路径）。添加你所在的文件夹

通过单击添加文件夹按钮保存了文件。选择文件夹，然后单击

按钮保存→关闭。

（Set the path to MATLAB to run the file (File → Set Path). Add the folder in which you
have saved the file by clicking the button Add Folder. Choose the folder and click the
button Save → Close.）

3) 将以下代码添加到M-File中：

（Add the following codes to the M-File:）

4) 保存文件。（文件→另存为）（Save the file. (File → Save As)

5）通过单击调试→运行filename.m运行M-File。（Run the M-File by clicking Debug → Run filename.m.）

6) 保存数字。（Save the figures.）

7) 将绘图命令更改为茎和楼梯，然后运行文件。（Change the plot command to stem and stairs and run the file.）

Alright,Let's do some tasks now😊

Task 1:

使用下面的公式创建葡萄糖摄入量(𝐺𝑛)图。复制该方程并进行一些修改，生成如图1所示的图。
(Use the equation below to create Glucose intake (𝐺𝑖𝑛) graph. Duplicate the equation with some modification to produce the graph as in Figure 1.)

 

           Figure 1

 

Task 1 Results Code:

% Define the time intervals t = linspace(0, 400, 4000);

% Initialize Gin
Gin = zeros(size(t));
% Loop over the time intervals

for i = 1:length(t)
     if t(i) >= 0 && t(i) < 17.5
         Gin(i) = 0.25 + (4.75/17.5)*t(i); % Increasing part of the first spike
     else if t(i) >= 17.5 && t(i) < 35
         Gin(i) = 5 - (4.75/17.5)*(t(i)-17.5); % Decreasing part of the first spike
     else if t(i) >= 35 && t(i) < 240
         Gin(i) = 0.25; % Constant part between the spikes
     else if t(i) >= 240 && t(i) < 257.5
         Gin(i) = 0.25 + (4.75/17.5)*(t(i)-240); % Increasing part of the second spike
     
     else if t(i) >= 257.5 && t(i) < 275
        Gin(i) = 5 - (4.75/17.5)*(t(i)-257.5); % Decreasing part of the second spike
     
     else
        Gin(i) = 0.25; % Constant part after the second spike
      end

end
% Plot the graph
figure(1);
plot(t, Gin);
title('Glucose Intake Rate');

xlabel('Time (t)');

ylabel('Gin(t)');
grid on;

 

                                     Task 1 Result 

 

Task 2:

使用下面的公式创建一个Lispro胰岛素摄入量(𝐼𝑛)图。复制该方程并进行一些修改，生成如图2所示的图。

(Use the equation below to create a Lispro insulin intake (𝐼𝑖𝑛) graph. Duplicate the equation with some modification to produce the graph as in Figure 2.)

Figure 2

 

Task 2 Result Code:

x_values = 0:1:480;
y = zeros(size(x_values));

pattern_length = 240;
for i = 1:length(x_values)
      x = mod(x_values(i), pattern_length); % Use modulo to repeat the pattern
if (x >= 0 && x < 15)

       y(i) = 0.25;

elseif (x >= 15 && x < 30)
        y(i) = 0.25 + 2*(1 + (x-30)/15);
elseif (x >= 30 && x < 120)
        y(i) = 0.25 + 2*(1 - (x-30)/90);
elseif (x >= 120 && x < pattern_length)

         y(i) = 0.25;
  end

end
figure;
plot(x_values, y);

title('Insulin intake rate');

xlabel('Time');

ylabel('Value'); 

 

Task 2 Result 

 

2.0 脚本-文件和绘图(Script - file and plotting)

1) 单击文件→新建→空白M文件打开一个新的M文件。(Open a new M-file by clicking File → New → Blank M-File.)

2) 将以下代码添加到M-File中: (Add the following codes to the M-File)

3) 保存文件。（文件→另存为）(Save the file. (File → Save As))

4) 通过单击调试→运行filename.m运行M-File。(Run the M-File by clicking Debug → Run filename.m.)

Let's us do related two tasks

 Task 1:

在MATLAB中绘制函数（请参阅图形绘制线条和标记）

● y=x^2 −2x–3,x = [-10:10]带钻石标记

● y = 3xe^−0.25，x = [−5:25] 带圆形标记

● y = 1/(1 + 99e^−0.5t)，t = [0:30] 带有洋红色线

[ Plot the functions in MATLAB (See Graphics plotting lines and marker)
● y = 𝑥^2 −2x–3, x = [-10:10]with diamond marker
● y = 3𝑥𝑒^−0.25,x = [−5:25] with circle marker
● y = 1/(1 + 99𝑒^−0.5𝑡), t = [0:30] with magenta line ]

Table 1 : Symbol of Colour, Line Style and Marker in MATLAB

Task 1 Result Code :

% 1. Plot the functions in MATLAB

x = -10:10;
y = x.^2 - 2.*x - 3;
subplot(3,1,1)
plot(x, y, 'kd')
xlabel('x = [-10:10]')
ylabel('y = x^2 - 2x - 3')
title('Plot of the function y = x^2 - 2x - 3','FontSize',12)
x = -5:25;
y = 3.*x.*exp(-0.25.*x);
subplot(3,1,2)
plot(x, y, 'ko')
xlabel('x = [-5:25]')
ylabel('y = 3xe^-0.25x')
title('Plot of the function y = 3xe^-0.25x','FontSize',12)
t = 0:30;
y = 1./(1 + 99.*exp(-0.5.*t));
subplot(3,1,3)
plot(t, y, 'm-')
xlabel('t = [0:30]')
ylabel('y = 1/(1 + 99e^-0.5t)')
title('Plot of the function y = 1/(1 + 99e^-0.5t)','FontSize',12)-

 

 Result 1

Task 2:

在单个窗口中绘制功能

Y1 = (x − 2)(3 − x)(x + 5) 和 y2 = −2(x−2)(3−x)(x+5) with x = [−10:10] y1 = (x + 5)(3 − x) 和 y2 =−3(x + 5)(3 − x) with x = [−10:10]

[ Plot the functions in a single window
y1 = (x − 2)(3 − x)(x + 5) and y2 = −2(x−2)(3−x)(x+5) with x = [−10:10] y1 = (x + 5)(3 − x) and y2 =−3(x + 5)(3 − x) with x = [−10:10] ]

Task 2 Result Code:

% 2. Plot the functions in a single window x = -10:10;
y1 = (x - 2).*(3 - x).*(x + 5);
y2 = -2.*(x - 2).*(3 - x).*(x + 5);

subplot(2,1,1)
plot(x, y1, 'r-', x, y2, 'b-')
xlabel('x = [-10:10]')
ylabel('y1 and y2')
legend('y1 = (x - 2)(3 - x)(x + 5)', 'y2 = -2(x - 2)(3 - x)(x + 5)') title('Plot of the functions y1 and y2','FontSize',12)
y1 = (x + 5).*(3 - x);
y2 = -3.*(x + 5).*(3 - x);
subplot(2,1,2)
plot(x, y1, 'r-', x, y2, 'b-')
xlabel('x = [-10:10]')
ylabel('y1 and y2')
legend('y1 = (x + 5)(3 - x)', 'y2 = -3(x + 5)(3 - x)') title('Plot of the functions y1 and y2','FontSize',12)

 

Result 2

3.0 功能-文件（M-文件）(Function - File (M-File))

1) 打开一个新的M文件，通过单击文件→新建→空白M文件来创建函数文件。（Open a new M-File to create a function-file by clicking File → New → Blank M-File.）

2）将以下代码添加到M-File并保存：（Add the following codes to the M-File and save:）

3）将文件保存为func1.m。（文件→另存为）（Save the file as func1.m. (File → Save As)）

4）通过单击文件→新建→空白M文件创建脚本文件，并添加以下内容代码：（Create a script-file by clicking File → New → Blank M-File and add the following
codes:）

5）将文件保存为mainprog.m。（文件→另存为）（Save the file as mainprog.m. (File → Save As)）

6）单击调试→运行mainprog.m运行脚本文件（Run the script-file by clicking Debug → Run mainprog.m）

 

Task 1

编写一个计算向量标准差的函数文件。该函数应该在向量x中读取，然后创建一个脚本文件来调用该函数并获得平均值。标准差方程：

其中n=length(x)，x(i)是x的第i个分量（我们数字列表中的第i个条目），xaverage是平均值。

创建您自己的矢量x值。使用脚本文件从主脚本调用标准偏差函数。

提示：在Matlab中使用长度、总和和sqrt函数。

不要在程序中使用内置的std()函数。然而，您的函数的结果应该与内置的std()函数相同。
 

Task 1 Result Code:

Function script
% std_deviation.m
function s = std_deviation(x)
        n = length(x);

        x_average = sum(x) / n;
        squared_diff = (x - x_average).^2;
        s = sqrt(sum(squared_diff) / (n - 1));

end
 

Main script
% Create your own vector x x = [2, 4, 4, 4, 5, 5, 7, 9];
% Call the custom standard deviation function custom_std_result = std_deviation(x);
% Call the built-in std() function built_in_std_result = std(x);
% Display the results
disp(['Standard Deviation (Custom Function): ', num2str(custom_std_result)]);

disp(['Standard Deviation (Built-in std() Function): ', num2str(built_in_std_result)]);

Result 1

 

Task 2

编写一个函数文件，将温度以华氏度（◦F）转换为摄氏度（◦C）为单位。使用输入命令让用户输入值和fprintf命令来显示文本和数字的混合，例如“温度为100摄氏度”。回想一下转换公式，C = 5/9 ∗ (F − 32)。

Task 2 Result Code:

Function script
function celsius = fahrenheit_to_celsius(fahrenheit)
    % Check if the user provided the Fahrenheit temperature

if nargin < 1
    error('Please provide a temperature in degrees Fahrenheit.');

end
% Convert temperature from Fahrenheit to Celsius
celsius = (5/9) * (fahrenheit - 32);
% Display the result using fprintf
fprintf('The temperature is %.2f degrees Celsius.\n', celsius);
end
 

Main script
% Example usage
temperature_F = input('Enter temperature in degrees Fahrenheit: '); fahrenheit_to_celsius(temperature_F);

 

 

4.0 信号生成和卷积

1)通过单击文件→新建→空白M文件打开新的M文件。(Open a new M-File by clicking File → New → Blank M-File.)

2)保存文件。（文件→另存为）(Save the file. (File → Save As))

3)设置MATLAB运行文件的路径（文件→设置路径）。单击“添加文件夹”按钮，添加您已保存文件的文件夹。选择文件夹，然后单击保存→关闭按钮。(Set the path for MATLAB to run the file (File → Set Path). Add the folder in which you have saved the file by clicking button Add Folder. Choose the folder and click the button Save → Close.)

4)将以下代码添加到M-File中：(Add the following codes to the M-File:)

 

5)单击调试→运行filename.m来运行M文件(Run the M-File by clicking Debug → Run filename.m)

Task 1:

 

为下图绘制带有输入x（t）和脉冲响应h（t）的输出y（t）：

A）输入x（t）具有矩形形式，脉冲响应h（t）具有正方形形式：

（使用持续时间为-1:0.01:5。提示：应用图（t，y（1：长度（t）））

 

Task 1 result Code:

%% Convolution Task 1 %% %%%%%%%%%%%%%%%%%%%%%%%%% tint = -1;
tfinal = 5;
tstep = 0.01;
t = tint:tstep:tfinal;
% a) Create a rectangle waveform from t = -1 to 0 x = 1 * ((t >= 0) & (t <= 2));
% Shift the rectangle waveform to the right (from 0 to 1) x_shifted = circshift(x, [0, length(t) - 1]);
% b) Create an impulse response h(t) with a square form from t = 0 to 1 h = 1 * ((t >= 0) & (t <= 1));
% c) Convolution of x(t) and h(t)
t2 = tint: tstep: tfinal; % shift the range from 1 to 4
y = conv(x_shifted, h) * tstep; % convolution using conv function is in discrete
% Plotting
figure(2) subplot(3,1,1)
plot(t, x)
axis([-1 5 0 2]) title('Input Signal x(t)')
subplot(3,1,2)
plot(t, h)
axis([-1 5 0 2])
title('Impulse Response h(t)')
subplot(3,1,3)
plot(t2, y(1:length(t2)))
axis([-1 5 0 2])
title('Convolution of x(t) and h(t) shifted to the right')

Task 1 Result 

 

Task 2:

输入x（t）具有三角形形式，脉冲响应h（t）具有矩形形式：

 

Task 2 Result Code: 

tint = 0;
tfinal = 5;
tstep = 0.01;
t = tint:tstep:tfinal;
% create a rectangle waveform from t = 0-1 x = t .* (t >= 0 & t <= 1);
figure(4)
subplot(3,1,1)
plot(t,x)
axis([0 5 0 1.5])
title('Triangle Waveform')
% create a rectangle waveform from t = 0-1 h = 1*((t>=0)&(t<=3));
subplot(3,1,2)
plot(t,h)
axis([0 5 0 1.5]) title('Rectangle Waveform')

t2=2*tint:tstep:2*tfinal; % output has double length of input % convolution of 2 waveform
y = conv(x,h)*tstep; % convolution using conv function is in discrete % therefore need to multiply with delta t/step size subplot(3,1,3),plot(t2,y)
axis([0 5 0 1.5])
title('Convolution of 2 Waveforms')

 

Task 2 Result

 

5.0 傅里叶变换

1) 通过单击文件→新建→空白M文件打开新的M文件。(Open a new M-File by clicking File → New → Blank M-File.)

2) 保存文件。（文件→另存为）(Save the file. (File → Save As))

3) 设置MATLAB运行文件的路径（文件→设置路径）。单击“添加文件夹”按钮，添加您保存文件的文件夹。选择文件夹，然后单击保存→关闭按钮。(Set the path for MATLAB to run the file (File → Set Path). Add the folder in which you have saved the file by clicking the button Add Folder. Choose the folder and click the button Save → Close.)

4) 将以下代码添加到M-File中：(Add the following codes to the M-File:)

 

5) 单击调试→运行filename.m来运行M文件(Run the M-File by clicking Debug → Run filename.m)

 

Task a and Task b

 Task a Result Code:

% Task Exp 5 % Part a)

fs = 8000;
N = 1000;
n = 0:N-1;
x = 2 * sin(2000 * pi * n / fs);

% Compute FFT
X = fft(x, N);
% Frequency axis
f = linspace(0, fs, N);
% Plot the magnitude spectrum
figure;
plot(f, abs(X));

title('Task 1');

xlabel('Frequency (Hz)');

ylabel('Magnitude');

Task a Result 

 

Task b Result Code:

% Task Exp 5
% Part b)
n = 0:29; % Time range for the second signal
% Values of N
N1 = 64;
N2 = 128;
N3 = 256;
N4 = 512;
N5 = 1024;
% Signal definition
x = cos(2 * pi * n / 10);
% Compute FFT for each value of N
X1 = fft(x, N1);
X2 = fft(x, N2);
X3 = fft(x, N3);
X4 = fft(x, N4);
X5 = fft(x, N5);
% Frequency axis
f1 = linspace(0, fs, N1);
f2 = linspace(0, fs, N2);
f3 = linspace(0, fs, N3);
f4 = linspace(0, fs, N4);
f5 = linspace(0, fs, N5);
% Plot the magnitude spectrum for different values of N figure;
subplot(3, 2, 1); plot(f1, abs(X1)); title('N = 64'); subplot(3, 2, 2); plot(f2, abs(X2)); title('N = 128'); subplot(3, 2, 3); plot(f3, abs(X3)); title('N = 256');

subplot(3, 2, 4); plot(f4, abs(X4)); title('N = 512'); subplot(3, 2, 5); plot(f5, abs(X5)); title('N = 1024');

Task b Result

6.0 滤波器的实现

设备：

a）MATLAB软件

b)  心电图垫文件/信号处理工具箱（ECG mat file / signal processing toolbox）

1) 将以下代码写入新的M文件：

2) 单击调试→运行filename.m来运行M文件 (Run the M-File by clicking Debug → Run filename.m)

Task 

继续下面的编码，对嘈杂的心电图（心电图）信号实施滤波器：(Continue the coding below with implementation of filter to the noisy ECG (Electrocardiogram) signal:

 Task Result Code:

Ecg code
function x = ecg(L)
a0=[0, 1,40, 1, 0,-34,118,-99, 0, 2, 21, 2, 0, 0, 0];
d0 = [0, 27, 59, 91, 131, 141, 163, 185, 195, 275, 307, 339, 357, 390, 440]; a = a0 / max(a0);
d = round(d0 * L / d0(15));
d(15) = L;
for i = 1:14
   m = d(i) : d(i+1) - 1;
   slope = (a(i+1) - a(i)) / (d(i+1) - d(i));
   x(m+1) = a(i) + slope * (m - d(i)); %#ok<AGROW>
end

Filter Main Code
% Load ecg and simulate noisy ECG Fs = 500;
x = repmat(ecg(Fs), 1, 8);
x = x + randn(1, length(x)).*0.18;
% Plot noisy signal figure
subplot(211), plot(x) title('Noisy ECG')
b = ones(1, 20)/20;
%Implementation of filter y = filter(b,1,x); subplot(212), plot(y) xlabel('Time') ylabel('Amplitude') title('Filtered ECG')

Task Result