Biosignal and Medical Image Processing, 3e

• Written for students and bioengineers, Biosignal and Medical Image Processing presents a range of signal and image processing tools for biomedical applications. The book provides a working knowledge of the technologies addressed while highlighting important implementation procedures, common pitfalls, and essential application concepts. Topics include spectral analysis, digital filters, wavelet transforms, multivariate analysis, and image segmentation and reconstruction.
• MATLAB and Image Processing Toolbox are used to solve application examples in the book.
• https://www.crcpress.com/Biosignal-and-Medical-Image-Processing-Third-Edition/Semmlow-Griffel/9781466567368
• http://www.mathworks.com/support/books/book91715.html?category=14&language=1&view=category
  

1. % Ex 2.1 Digital solution. Generate a 1 cycle sine wave.  Assume a sample interval
   
2. %  of Ts = 0.005 sec. and make N = 500 points long (both arbitrary).
   
3. %  So solving Eq. 1.4 for the total time: TT = NTs = 0.005(500) = 2.5 sec.
   
4. %  To generate cycle given these parameters, the frequency of the sine wave  
   
5. %  should be f = 1/TT = 1/2/5 = 0.4 Hz.   
   
6. %
   
7. close all; clear all;
   
8. N = 500;                 % Number of points for waveform
   
9. Ts = .005;                                 % Sample interval =  5 msec
   
10. t  = (1:N)*Ts;                     % Generate time vector (t = N Ts)
    
11. f = 1/(Ts*N);            % Sine wave freq for 1 cycle
    
12. A = 1;                   % Sine wave amplitude
    
13. x = A*sin(2*pi*f*t);       % Generate sine wave
    
14. plot(t,x);               % Plot (not shown)
    
15. rms = sqrt(mean(x.^2))   % Evaluate the RMS value and output

复制代码

1. % Example 2.2 Evaluate a signal to determine if it is stationarity.
   
2. % If not, attempt to modify the signal to remove the nonstationarity if
   
3. % possible.  The file data_c1.mat contains the signal in variable x.  
   
4. % The signal is 1000 points long and the sample interval is Ts = 0.001 sec.  
   
5. %
   
6. %Solution, Part 1: In Part 1 of this example we want to determine this signal
   
7. % is stationary.  If it is not, we will modify the signal in Part 2 to make
   
8. % it approximately stationary. There are a number of approaches we could try,
   
9. % but one straightforward method would be to segment the data and evaluate the
   
10. % mean and variance of the segments.  If they are the same for all segments,
    
11. % the data is probably stationary.  If these measures change segment-to-segment,
    
12. % the signal is clearly nonstationary.
    
13. close all; clear all;
    
14. % Example 2.2, Part 1. Evaluate a signal for stationarity.
    
15. load data_c1;                   % Load the data. Signal in variable x
    
16. for k = 1:4                     % Segment signal into 4 segments
    
17.     m = 250*(k-1) + 1;          % Index of first segment sample
    
18.     segment = x(m:m+249);       % Extract segment
    
19.     avg(k) = mean(segment);     % Evaluate segment mean
    
20.     variance(k) = var(segment); %   and segment variance
    
21. end
    
22. disp('Mean   Segment 1 Segment 2 Segment 3 Segment 4')  % Display heading
    
23. disp(avg)                                                                  % Output means
    
24. disp('Variance Segment 1 Segment 2 Segment 3 Segment 4')% Heading
    
25. disp(variance)                                                          % Output variance       
    
26. %
    
27. % Solution, Part 2:  One trick is to transform the signal by taking the
    
28. % derivative: just taking the difference between each point.  
    
29. % Example 2.8, Part 2. Modify a nonstationary signal to become stationary
    
30. %
    
31. y = [diff(x);0];                % Take difference between points.
    
32. for k = 1:4                     % Segment signal into 4 segments
    
33.     m = 250*(k-1) + 1;          % Index of first segement sample
    
34.     segment = y(m:m+249);       % Extract segment
    
35.     avg(k) = mean(segment);     % Evaluate segment mean
    
36.     variance(k) = var(segment); %   and segment variance
    
37. end
    
38. disp('Mean   Segment 1 Segment 2 Segment 3 Segment 4')  % Display heading
    
39. disp(avg)                                                                  % Output means
    
40. disp('Variance Segment 1 Segment 2 Segment 3 Segment 4')% Heading
    
41. disp(variance)                                                          % Output variance       

复制代码

1. % Example 2.3 and Figure 2.4 Evaluation of the distribution of data produced by MATLAB's
   
2. % rand and randn functions.  
   
3. %
   
4. clear all; close all;
   
5. N = 20000;                          % Number of data points
   
6. nu_bins = 40;                   % Number of bins
   
7. y = randn(1,N);                     % Generate random Gaussian noise
   
8. [ht,xout] = hist(y,nu_bins);         % Calculate histogram
   
9. ht = ht/max(ht);                % Normalize histogram to 1.0
   
10. subplot(2,1,1);
    
11. bar(xout, ht);                     % Plot as bar graph
    
12. xlabel('x','FontSize',14);
    
13. ylabel('P(x)','FontSize',14);
    
14. title('Gaussian Distribution','Fontsize',14);
    
15. % Repeat for uniform distribution
    
16. y = rand(1,N);                     % Generate random Gaussian noise
    
17. [ht,xout] = hist(y,nu_bins);         % Calculate histogram
    
18. ht = ht/max(ht);                % Normalize histogram to 1.0
    
19. subplot(2,1,2);
    
20. bar(xout, ht);                     % Plot as bar graph
    
21. xlabel('x','FontSize',14);
    
22. ylabel('P(x)','FontSize',14);
    
23. title('Uniform distribution','Fontsize',14);

复制代码

1. % Example 2.4 This example evaluates the noise reduction provided by ensemble
   
2. % averaging. An ensemble of simulated visual evoked responses is found as
   
3. % matrix variable ver in MATLAB file ver.mat.  The actual, noise-free VER
   
4. % signal is found as variable actual_ver in the same file.  Ts = 0.005 sec.
   
5. % Construct an ensemble average of 100 responses.  Subtract out the noise-free
   
6. % VER from a selected individual VER and form the ensemble average to get
   
7. % the noise in each record.  Compare the reduction in standard deviation
   
8. % produced by averaging with that predicted theoretically.
   
9. %
   
10. load ver;                       % Get visual evoked response data;
    
11. Ts = 0.005;                     % Sample interval = 5 msec
    
12. [nu,N] = size(ver);             % Get data matrix size
    
13. if nu > N
    
14.     ver = ver';
    
15.     t = (1:nu)*Ts;              % Generate time vector
    
16. else
    
17.     t = (1:N)*Ts;               % Time vector if no transpose   
    
18. end
    
19. %
    
20. subplot(2,1,1);
    
21. plot(t,ver(4,:),'k');                     % Plot individual record (Number 4)
    
22.   xlabel('Time (sec)','FontSize',14);
    
23. ylabel('VER','FontSize',14);
    
24. % Construct and plot the ensemble average
    
25. avg = mean(ver);                % Take ensemble average
    
26. subplot(2,1,2);
    
27. plot(t,avg,'k');                          % Plot ensemble average other data
    
28. xlabel('Time (sec)','FontSize',14);
    
29. ylabel('VER','FontSize',14);
    
30. %
    
31. % Estimate the noise components
    
32. unaveraged_noise = ver(4,:) - actual_ver; % Noise in signal
    
33. averaged_noise = avg - actual_ver;        % Noise in ensemble avg.
    
34. std_unaveraged = std(unaveraged_noise);   % Std or noise in signal
    
35. std_averaged = std(averaged_noise);             % Std of noise in ensemble avg.
    
36. theoretical =  std_unaveraged/sqrt(100);  % Theoretical noise reduction
    
37. disp(' Unaveraged Averaged Theroetical')        % Output results
    
38. disp([std_unaveraged std_averaged theoretical])

复制代码

1. % Example 2.5 Find the angle between two vectors in deg: x = [1.7, 3, 2.2] and
   
2. % y = [2.6, 2.4, 3.2]. Since these vectors are short, plot the two in 3-Dimensions.
   
3. % Solution.  Consruct the two vectors and find the scalar product using
   
4. % Eq. 2.32 (same as Eq.2.31).  Solve Eq. 2.33 for the angle and apply.
   
5. %
   
6. clear all; close all;
   
7. x = [1.7, 3, 2.2];              % Generate the vectors
   
8. y = [2.6, 1.6, 3.2];
   
9. sp = sum(x.*y);                 % Calculate the scalar product
   
10. mag_x = sqrt(sum(x.^2));        % Calculate magnitude of x vector
    
11. mag_y = sqrt(sum(y.^2));        % Calculate magnitude of y vector
    
12. cos_theta = sp/(mag_x*mag_y);   % Calculate the cosine of theta
    
13. angle = acos(cos_theta);        % Take the arc cosine and
    
14. angle = angle*360/(2*pi);       %    convert to degrees
    
15. disp(angle)                     % Display result
    
16. plot3(x(1),x(2),x(3),'k*','LineWidth',2);      % Plot x vector end point
    
17. hold on;
    
18. plot3([0 x(1)],[0 x(2)],[0 x(3)],'k','LineWidth',2);   % Plot x vector line
    
19. plot3(y(1),y(2),y(3),'k*','LineWidth',2);      % Plot y vector end point
    
20. plot3([0 y(1)],[0 y(2)],[0 y(3)],'k','LineWidth',2);   % Plot vector y line
    
21. title(['Angle = ',num2str(angle,2),' (deg)'],'FontSize',14);  % Output angle
    
22. grid on;
    
23. xlabel('x','FontSize',16); ylabel('y','FontSize',16); zlabel('z','FontSize',16);

复制代码

1. % Example 2.6 Generate a 500 point, 2 Hz sine wave and a 4
   
2. % Hz sine wave of the same length.  Make TT = 1 sec. Are these two waveforms orthogonal?
   
3. % Solution. Generate the waveforms. Since N = 500 and  TT = 1 sec,
   
4. % Ts = TT /N = 0.002 sec.  Apply Eq. 2.29 to these to waveforms.  
   
5. % If the result is small (near zero) the waveforms are orthogonal.
   
6. %
   
7. clear all close all;
   
8. Ts = 0.002;             % Sample interval
   
9. N = 500;                % Number of points
   
10. t = (0:N-1)*Ts;         % Time vector
    
11. f1 = 2;                 % Frequency of sine wave 1
    
12. f2 = 4;                 % Frequency of sine wave 2
    
13. x = sin(2*pi*f1*t);     % Sine wave 1
    
14. y = sin(2*pi*f2*t);     % Sine wave 2
    
15. r = sum(x.*y);          % Eq. 2.29
    
16. disp(r)

复制代码