11.1 - BME 240 - Spring2024

pdf

School

University of Illinois, Chicago *

*We aren’t endorsed by this school

Course

240

Subject

Electrical Engineering

Date

Apr 3, 2024

Type

pdf

Pages

49

Uploaded by BaronBraverySquid16

Report
Signal analysis in the frequency domain Week 11.1 - BME 240 Spring2024
BB: TAEMG3.txt data Sample frequency was 500 Hz 1. Plot the frequency spectrum of raw EMG signals 2. Remove the frequency at 60 Hz (58-62Hz) using band stop filter (Butter, 4 th order); 3. Remove the frequency at 180 Hz (178-182 Hz) using band step (Butter 4 th order). 4. Plot the frequency spectrum after band stop filter. 5. Plot the EMG after removing the noise.
[newfile, newpath] = uigetfile('*.*', 'Select result File'); cd(newpath); filenameH2='TAemg3.txt'; dataH_2=load(filenameH2); figure(1); figure;plot(dataH_2); N1=length(dataH_1); fe1=fft(dataH_1); f1=(1:N1)*fs/N1; figure;plot(f1(1:(N1/2)-1),abs(fe1(2:(N1/2))));
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
0 50 100 150 200 250 0 50 100 150 200 250 0 50 100 150 200 250 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 X 180 Y 1.91192 What kind of filter should we use?
filtfilt Zero-phase forward and reverse digital IIR filtering. Y = filtfilt (B, A, X) filters the data in vector X with the filter described by vectors A and B to create the filtered data Y. filter One- dimensional digital filter. Y = filter (B,A,X) filters the data in vector X with the filter described by vectors A and B to create the filtered data Y.
butter Butterworth digital and analog filter design. [B, A] = butter (N, Wn) designs an Nth order lowpass digital Butterworth filter and returns the filter coefficients in length N+1 vectors B (numerator) and A (denominator). The cutoff frequency Wn must be 0.0 < Wn < 1.0, with 1.0 corresponding to half the sample rate fs/2. [B,A] = butter(N,Wn,' high ') designs a highpass filter. [B,A] = butter(N,Wn,' low ’) designs a lowpass filter. [B,A] = butter(N,Wn,' stop ') is a bandstop filter if Wn = [W1 W2].
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
WN=[55/(fs/2) 65/(fs/2)]; [b,a]=butter(4,WN, 'stop' ); xemgn=filtfilt(b,a,xemg); figure;plot(xemgn); 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 -0.15 -0.1 -0.05 0 0.05 0.1 0.15 0 500 1000 1500 2000 2500 3000 3500 4000 4500 5000 -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 Band stop filter
xemg=xemgn; N=length(xemg); xfe=fft(xemg); fs=500; T=N/fs; f=(1:N)*fs/N; figure;plot(f(1:(N/2)-1),abs(xfe(2:(N/2)))) 0 50 100 150 200 250 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
wn=[55/(fs/2) 65/(fs/2)]; [b,a]=butter(4,wn, 'stop' ); xemgn1=filtfilt(b,a,dataH_1(:)); figure;plot(xemgn1); fe2=fft(xemgn1); figure;plot(f1(1:(N1/2)-1),abs(fe2(2:(N1/2)))); wn=[178/(fs/2) 182/(fs/2)]; [b,a]=butter(4,wn, 'stop' ); xemgn2=filtfilt(b,a,xemgn1); fe3=fft(xemgn2); figure;plot(f1(1:(N1/2)-1),abs(fe3(2:(N1/2)))); figure;plot(xemgn2);
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
wn=[10/(fs/2)]; [b,a]=butter(4,wn, 'high' ); xemgn3=filtfilt(b,a,xemgn2); figure;plot(xemgn3); fe5=fft(xemgn3); figure;plot(f1(1:(N1/2)-1),abs(fe5(2:(N1/2)))); Remove low frequency noise What kind of filter should we use?
EMG after high-pass filtering with cutoff frequency at 10Hz
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
wn=[20/(fs/2)]; [b,a]=butter(4,wn, 'high' ); xemgn3=filtfilt(b,a,xemgn2); figure;plot(xemgn3); fe5=fft(xemgn3); figure;plot(f1(1:(N1/2)-1),abs(fe5(2:(N1/2))));
Example EMG signals with electrical stimulation induced artifacts
[newfile, newpath] = uigetfile( '*.*' , 'Select result File' ); cd(newpath); filenameH5= 'T11_80Hz' ; dataH_5=load(filenameH5); figure; subplot(3,1,1);plot(dataH5(5,:)); subplot(3,1,2);plot(dataH5(6,:)); subplot(3,1,3);plot(dataH5(7,:));
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
xemg=dataH(4,:); N=length(xemg); xfe=fft(xemg); fs=500; T=N/fs; f=(1:N-1)*fs/N; figure;plot(f(1:(N/2)-1),abs(xfe(2:(N/2))));
WN=[75/(fs/2) 85/(fs/2)]; [b,a]=butter(4,WN, 'stop' ); dataH1(4,:)=filtfilt(b,a,dataH(4,:)); WN=[235/(fs/2) 245/(fs/2)]; [b,a]=butter(4,WN, 'stop' ); dataH2(4,:)=filtfilt(b,a,dataH1(4,:)); WN=[95/(fs/2) 105/(fs/2)]; [b,a]=butter(4,WN, 'stop' ); dataH3(4,:)=filtfilt(b,a,dataH2(4,:)); WN=[55/(fs/2) 65/(fs/2)]; [b,a]=butter(4,WN, 'stop' ); dataH4(4,:)=filtfilt(b,a,dataH3(4,:));
[b1,a1]=butter(4,20/(fs/2), 'high' ); % low pass at 20 Hz dataH5(4,:)=filtfilt(b1,a1,dataH4(4,:)); After removed the noises
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
[newfile, newpath] = uigetfile( '*.*' , 'Select result File' ); cd(newpath); filenameH5= 'position.txt' ; dataH_5=load(filenameH5); figure;plot(dataH_5);
wn=10/250; [b1,a1]=butter(4,wn, 'low' ); xpe=filtfilt(b1,a1,dataH_1); xpen=fft(xpe); figure;plot(f(1:(N/2)-1),abs(xpen(2:(N/2)))) figure;plot(xpe); Low pass filter
xpe=filtfilt(b1,a1,dataH_5); xpen=fft(xpe); figure;plot(fp(1:(Np/2)-1),abs(xpen(2:(Np/2)))) figure;plot(xpe); 0 50 100 150 200 250 0 5000 10000 15000 0 1 2 3 4 5 6 7 8 9 10 0 1000 2000 3000 4000 5000 6000 7000 8000
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
filenameH6= 'vertforce.txt' ; dataH_6=load(filenameH6); figure;plot(dataH_6);
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
force=dataH_6; N=length(force); xfe=fft(force); fs=500; T=N/fs; f=(1:N-1)*fs/N; figure;plot(f(1:(N/2)-1),abs(xfe(2:(N/2))));
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
[b1,a1]=butter(4,25/(fs/2), 'low' ); % low pass at 25 Hz forcen=filtfilt(b1,a1,force); figure;plot(forcen);
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
dxpe=diff(xpe,1)/(1/fs); figure;plot(dxpe); diff Difference and approximate derivative. diff (X), for a vector X, is [ X(2)-X(1) X(3)-X(2) ... X(n)-X(n-1)].
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
figure;plot(forcen); df=diff(forcen,1); figure;plot(df);
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
emg1=interp(signal(peakind(1):peakind(2),1),100); ln=length(signal(peakind(1):peakind(2))) length(emg1) figure;plot(emg1); emg2=downsample(emg1,ln); length(emg2)
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
Synchronized position signals
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
[b1,a1]=butter(4,10/(fs/2), 'low' ); % low pass at 10 Hz p1=filtfilt(b1,a1,position);
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
step=800; shift=300; st=1; for i=1:60 [maxv,indv]=max(p1(st:st+step)); peak(i)=maxv; peakind(i)=indv+st; st=peakind(i)+shift; end figure;plot(p1); hold on ;plot(peakind,peak, 'r*' ); Find the peaks
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
xemgn1=dataH5(6,:); for j=1:10 emg1=interp(xemgn1(peakind(j):peakind(j+1)),300); lg=peakind(j+1)-peakind(j); emg1_d(j,:)=downsample(emg1,lg); end figure; subplot(1,1,1);plot(mean((emg1_d)));xlim([1 300]) Interp and downsample EMG signals
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
Averaged EMG signals (baseline)
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
si=15; for j=1:10 emg1=interp(xemgn1(peakind(j+si):peakind(j+1+si)),300); lg=peakind(j+1+si)-peakind(j+si); emg1_dsi(j,:)=downsample(emg1,lg); end figure; subplot(1,1,1);plot(mean((emg1_dsi)));xlim([1 300])
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
figure;plot(abs(mean(emg1_d))); hold on ;plot(-abs(mean(emg1_dsi)), 'r' ); Compare the differences emgcomp1=[sum(abs(mean(emg1_d))) sum(abs(mean(emg1_dsi)))] emgcomp1 = 0.2769 0.3003 emgcomp=[rms(abs(mean(emg1_d))) rms(abs(mean(emg1_dsi)))] area RMS emgcomp = 0.0016 0.0019
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
ankle_xyz_baseline.txt Load file xRa, yRa, zRA, thet2R Fs=500 Hz. Run fft for the first 3 column data Run low pass filter with cut off at 10 hz Plot the frequency spectrum of the position signals to confirm. Based on column 4 data to segment the data (max peaks, 31 steps) Interp 100 points for each cycle Downsample to 100 points for each cycle Average x, y, z across 30 cycles (multiply -1 to the average y). Plot3 ensemble-averaged profiles of these position data.
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
-200 -100 0 100 200 300 400 500 500 520 540 560 580 600 620 640 660 680 700
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
[newfile, newpath] = uigetfile( '*.*' , 'Select result File' ); cd(newpath); filenameH1= 'ankle_xyz_baseline.txt' ; dataH_1=load(filenameH1); fs=500; N=length(dataH_1(:,1)); figure; subplot(4,1,1);plot(dataH_1(:,1)); subplot(4,1,2);plot(dataH_1(:,2)); subplot(4,1,3);plot(dataH_1(:,3)); subplot(4,1,4);plot(dataH_1(:,4));
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
0 0.5 1 1.5 2 2.5 3 10 4 0 200 400 0 0.5 1 1.5 2 2.5 3 10 4 -700 -600 -500 0 0.5 1 1.5 2 2.5 3 10 4 300 350 400 0 0.5 1 1.5 2 2.5 3 10 4 50 100
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
fs=500; t=(1:N)/fs; f=(1:N)*fs/N; X1=fft(dataH_1(:,1)); X2=fft(dataH_1(:,2)); X3=fft(dataH_1(:,3)); figure;subplot(3,1,1);plot(f(1:N/2-1),abs(X1( 2 :N/2))); subplot(3,1,2);plot(f(1:N/2-1),abs(X2( 2 :N/2))); subplot(3,1,3);plot(f(1:N/2-1),abs(X3( 2 :N/2))); fft to get frequency spectrum
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
0 50 100 150 200 250 0 1 2 3 10 6 0 50 100 150 200 250 0 5 10 10 5 0 50 100 150 200 250 0 1 2 10 5 m=20*((N/2)/250) figure;subplot(3,1,1);plot(f(1:m-1),abs(X1(2:m))); subplot(3,1,2);plot(f(1:m-1),abs(X2(2:m))); subplot(3,1,3);plot(f(1:m-1),abs(X3(2:m))); 0 2 4 6 8 10 12 14 16 18 20 0 1 2 3 10 6 0 2 4 6 8 10 12 14 16 18 20 0 5 10 10 5 0 2 4 6 8 10 12 14 16 18 20 0 1 2 10 5
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
[b,a]=butter(4,10/(fs/2), 'low' ); dataH_1n(:,1)=filtfilt(b,a,dataH_1(:,1)); dataH_1n(:,2)=filtfilt(b,a,dataH_1(:,2)); dataH_1n(:,3)=filtfilt(b,a,dataH_1(:,3)); X1n=fft(dataH_1n(:,1)); X2n=fft(dataH_1n(:,2)); X3n=fft(dataH_1n(:,3)); figure;subplot(3,1,1);plot(f(1:m-1),abs(X1n(2:m))); subplot(3,1,2);plot(f(1:m-1),abs(X2n(2:m))); subplot(3,1,3);plot(f(1:m-1),abs(X3n(2:m)));
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
0 2 4 6 8 10 12 14 16 18 20 0 1 2 3 10 6 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 10 5 0 2 4 6 8 10 12 14 16 18 20 0 0.5 1 1.5 2 10 5 base frequency is 0.76 Hz Cadence =0.76*60 =46 steps/minute
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
signal=dataH_1(:,4); step=800; shift=300; st=1; for i=1:43 [maxv,indv]=max(signal(st:st+step)); peak(i)=maxv; peakind(i)=indv+st; st=peakind(i)+shift; end figure;plot(signal(:,:)); hold on ;plot(peakind,peak, 'r*' );
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
0 0.5 1 1.5 2 2.5 3 10 4 40 50 60 70 80 90 100 110
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
for j=1:30 x1d=interp(dataH_1n(peakind(j):peakind(j+1),1),100); y1d=interp(dataH_1n(peakind(j):peakind(j+1),2),100); z1d=interp(dataH_1n(peakind(j):peakind(j+1),3),100); lg=peakind(j+1)-peakind(j); x1dn(j,:)=downsample(x1d,lg); y1dn(j,:)=downsample(y1d,lg); z1dn(j,:)=downsample(z1d,lg); end x1dnave=mean(x1dn(:,:)); y1dnave=mean(y1dn(:,:)); z1dnave=mean(z1dn(:,:)); figure;plot(x1dnave,-y1dnave) figure;plot3(x1dnave,-y1dnave,z1dnave)
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
-200 -100 0 100 200 300 400 500 520 540 560 580 600 620 640 660 680 700
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
x1dnstd=std(x1dn(:,:)); y1dnstd=std(-y1dn(:,:)); figure;plot(x1dnave,-y1dnave) hold on ;plot(x1dnave+x1dnstd,-y1dnave-y1dnstd, '-.' ) hold on ;plot(x1dnave-x1dnstd,-y1dnave+y1dnstd, '-.' ) -200 -100 0 100 200 300 400 500 520 540 560 580 600 620 640 660 680 700
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help
Load file emgbase3.txt First column is TA emg, second column is MG emg, and third column is position signal. Sample frequency was 500 Hz. 1) Run fft to plot the frequency spectrum of TA and MG. 2) Run high pass 4 th order butter filtering to remove low frequency noises (cutoff 10 Hz) 3) Run band stop butter filter to remove noise 60 Hz. 4) Run fft to plot the frequency spectrum of filtered TA and MG. 5) For the third column data run low pass 4 th order butter filtering to remove high frequency noises (cutoff frequency 10 Hz). 6) Find the max peaks of filtered position signals (30 cycles, third column data ). 7) Segment the TA and MG based on the position data to different cycles. Run interp to insert 300 points of the TA and MG for each cycle. 8) Run downsample to reduce the length of the TA and MG of each cycle to 300 points. 9) Average the TA and MG signals across 30 cycles. 10) Plot the averaged TA and MG signals ICP
Your preview ends here
Eager to read complete document? Join bartleby learn and gain access to the full version
  • Access to all documents
  • Unlimited textbook solutions
  • 24/7 expert homework help