Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- function plot_Fourier1(inputSignal,fs,wlen)
- % fs = sampling frequency
- % inputSignal = single column vector of numerical data
- % wlen = window length (recomended to be an integer power of 2)
- % shorter windows give higher resolution and are less affected by
- % high-frequency noise, they do however tend to loose the windowing effect
- % if they are chosen too small wrt the length of inputSignal
- % To run the function on the desired EMG data, run this from the command
- % line:
- % >> plot_Fourier1(Kanaal1mVolt,1000,512); (don't forget to import
- % Kanaal1mVolt)
- signal = inputSignal;
- x = signal(1: (length(signal)-1)); % exclude last value (NaN) because signal processor can't deal with that...
- time = 1 : (length(signal)-1);
- time = time.';
- time = time/1000; %divide by sampling frequency
- %%%%%%%%%%%%%%%%%%%%%%Plot Raw Data%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %define analysis parameters for stft
- xlen = length(x); % length of the signal
- hop = wlen/4; % hop size (recomended to be power of 2)
- nfft = 2^nextpow2(xlen);
- figure(1);
- plot(time,x);
- set(gca, 'FontName', 'Arial', 'FontSize', 14);
- xlabel('Time, s');
- ylabel('Amplitude, mV');
- title('Time Domain Signal - Raw EMG');
- % apply short-time Fourier transform for spectral analysis
- % number of fft points (recomended to be power of 2)
- % perform STFT
- [S, f, t] = stft(x, wlen, hop, nfft, fs);
- K = sum(hamming(wlen, 'periodic'))/wlen;
- % take the amplitude of fft(x) and scale it, so not to be a
- % function of the length of the window and its coherent amplification
- S = abs(S)/wlen/K;
- if rem(nfft, 2)
- S(2:end, :) = S(2:end, :).*2;
- else
- S(2:end-1, :) = S(2:end-1, :).*2;
- end
- % convert amplitude spectrum to dB (min = -120 dB)
- S = 20*log10(S + 1e-6);
- % plot the spectrogram
- figure(2);
- surf(t, f, S);
- colormap(jet);
- shading interp;
- axis tight;
- box on;
- view(0, 90);
- handl = colorbar;
- set(handl, 'FontName', 'Arial', 'FontSize', 14,'Colormap',jet);
- xlabel('Time, s');
- ylabel('Frequency, Hz');
- title('Spectogram of Raw EMG');
- ylabel(handl, 'Magnitude, dB');
- %%%%%%%%%%%%%%%%%%Filter the Data%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
- %notch filter to remove noise from 50 Hz
- % Design a filter with a Q-factor of Q=35 to remove a 50 Hz tone from
- % system running at 1000 Hz.
- Wo = 50/(fs/2); BW = Wo/35;
- [b,a] = iirnotch(Wo,BW);
- x = filter(b,a,x); %use filtfilt to avoid phase shifts?
- % Get rid of the DC component (hp filter at 350 Hz)
- fc=350;
- [b,a] = butter(6,2*fc/fs,'high');
- x = filter(b,a,x);
- %filter out noise from the electrical setup, including higher harmonics
- %of the AC supply (aka evrthg over 450 Hz)
- fc=450;
- [b,a] = butter(6,2*fc/fs,'low');
- x = filter(b,a,x);
- % number of fft points (recomended to be power of 2)
- % perform STFT
- [S1, f, t] = stft(x, wlen, hop, nfft, fs);
- % define the coherent amplification of the window
- %'periodic' — This option is useful for spectral analysis because it
- % enables a windowed signal to have the perfect periodic extension implicit
- % in the discrete Fourier transform. When 'periodic' is specified, hamming
- % computes a window of length L + 1 and returns the first L points.
- K = sum(hamming(wlen, 'periodic'))/wlen;
- % take the amplitude of fft(x) and scale it, so not to be a
- % function of the length of the window and its coherent amplification
- S1 = abs(S1)/wlen/K;
- if rem(nfft, 2)
- S1(2:end, :) = S1(2:end, :).*2;
- else
- S1(2:end-1, :) = S1(2:end-1, :).*2;
- end
- % convert amplitude spectrum to dB (min = -120 dB)
- S1 = 20*log10(S1 + 1e-6);
- % plot the spectrogram
- figure(3);
- surf(t, f, S1);
- colormap(jet);
- shading interp;
- axis tight;
- box on;
- view(0, 90);
- handl = colorbar;
- set(handl, 'FontName', 'Arial', 'FontSize', 14, 'Colormap', jet);
- xlabel('Time, s');
- ylabel('Frequency, Hz');
- title('Spectrogram of Filtered EMG');
- ylabel(handl, 'Magnitude, dB');
- set(gca,'Ylim',[0 550]);
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement