How do I segment with stockwell transfrom (s transform) ?
7 次查看(过去 30 天)
显示 更早的评论
I have the mathematical formula of stockwell transform, but I don't know how to write it as matlab code. Fs=2000Hz
2 个评论
Bora Eryilmaz
2022-12-14
Depends on what kind of noise you are trying to get rid off. Typically, you can use a low-pass filter: https://www.mathworks.com/help/signal/ref/lowpass.html
采纳的回答
Mathieu NOE
2022-12-15
hello
this can be a starter... load signal, bandpass filter, look at spectral content and tune your filter to get best results
the sampling rate is not yet defined so I took by default Fs = 2000 Hz , which gave me around one heart pulse every second (60 BPM)
clc
clearvars
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filename = 'matlab.mat';
data = load(filename);
signal = (data.data);
Fs = 2000; % sampling rate is TBD Hz
dt = 1/Fs;
[samples,channels] = size(signal);
% time vector
time = (0:samples-1)*dt;
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
for ck = 1:channels
newsignal(:,ck) = decimate(signal(:,ck),decim);
Fs = Fs/decim;
end
signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;
%% band pass filter section %%%%%%
f_low = 30;
f_high = 80;
N_bpf = 2; % filter order
[b,a] = butter(N_bpf,2/Fs*[f_low f_high]);
signal_filtered = filtfilt(b,a,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = samples; % for maximal resolution take NFFT = signal length.
OVERLAP = 0.75; % overlap will be used only if NFFT is shorter than signal length
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
plot(time,signal,'b',time,signal_filtered,'r');grid on
title(['Time plot / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend('raw','filtered');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
[freq, spectrum_filtered] = myfft_peak(signal_filtered,Fs,NFFT,OVERLAP);
figure(2),plot(freq,20*log10(spectrum_raw),'b',freq,20*log10(spectrum_filtered),'r');grid on
df = freq(2)-freq(1); % frequency resolution
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
legend('raw','filtered');
function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).
% Linear averaging
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer percentage of overlap % (between 0 and 0.95)
[samples,channels] = size(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
s_tmp = zeros(nfft,channels);
s_tmp((1:samples),:) = signal;
signal = s_tmp;
samples = nfft;
end
% window : hanning
window = hanning(nfft);
window = window(:);
% compute fft with overlap
offset = fix((1-Overlap)*nfft);
spectnum = 1+ fix((samples-nfft)/offset); % Number of windows
% % for info is equivalent to :
% noverlap = Overlap*nfft;
% spectnum = fix((samples-noverlap)/(nfft-noverlap)); % Number of windows
% main loop
fft_spectrum = 0;
for i=1:spectnum
start = (i-1)*offset;
sw = signal((1+start):(start+nfft),:).*(window*ones(1,channels));
fft_spectrum = fft_spectrum + (abs(fft(sw))*4/nfft); % X=fft(x.*hanning(N))*4/N; % hanning only
end
fft_spectrum = fft_spectrum/spectnum; % to do linear averaging scaling
% one sidded fft spectrum % Select first half
if rem(nfft,2) % nfft odd
select = (1:(nfft+1)/2)';
else
select = (1:nfft/2+1)';
end
fft_spectrum = fft_spectrum(select,:);
freq_vector = (select - 1)*Fs/nfft;
end
4 个评论
Mathieu NOE
2022-12-16
yes , this was posted in the "review" section
Works great once, but you may need to transpose the input.
The output is complex, so to get the magnitude plot in the picture you'll need to take the mag() or abs() of the result.
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!