How do I segment with stockwell transfrom (s transform) ?

I have the mathematical formula of stockwell transform, but I don't know how to write it as matlab code. Fs=2000Hz

[Bora Eryilmaz]

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

[Cey]

Thank you very much for your help. @Bora Eryilmaz

Answer 1

[Mathieu NOE]

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

[Cey]

First of all thank you very much for your help. After filtering, I need to segment with stockwell transform, but I couldn't find the code in matlab. I forgot I was going to change my question for segmentation.

[Mathieu NOE]

there are at least these two Fex submissions for you :

Stockwell Transform (S-Transform) - File Exchange - MATLAB Central (mathworks.com)

Time frequency distribution of a signal using S-transform (stockwell transform) - File Exchange - MATLAB Central (mathworks.com)

[Cey]

I looked at the 1st code but because it says 1xN size I inversed my data and the result was complex.

[Mathieu NOE]

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.