Clear Filters
Clear Filters

Problem in coding of STFT manually please correct me

1 view (last 30 days)
Hello guys ,I have a small problem in coding of short time fourier transform
actually what i wanted to do is to calculate the fft of a windowed signal and later shift the window by one segment i.e.., for ex: if i calculate the fft of y = x(1:98).*hamming(98) and then i wanted to subtract one segment from the signal that is Y(1) from the already calculated fft of y and then add a segment that is Y(99) so now i get the windowed fft from x(1:99).*hamming(98). even though i know the manual implementation of STFT ,I wanted to try this to check if it is computationally of low cost.THE main problem starts now ,how can i subtract a Y(1), can anyone of you give me the mathematical expression of it ??? I have also attached my code which has been done so far.
clear all
close all
clc
fs = 5000;
t = (0:100)/fs;
omega1 = 2*pi*500;
omega2 = 2*pi*200;
w = hamming(90)';
nw = length(w);%window with length of 98
x = sin(omega1*t)+ cos(omega2*t);%signal
N = length(x);
freq = (0:N-1)*fs/N;%frequency axis
w_total = [];
output = zeros(size(x));
%%%calculation of stft
for k = 0:nw-1
s = 0;
for n = 0 :nw-1
s = s + x(n+1).*w(n+1).*exp(-2i * pi * n * k/ N);
end
output(k+1) = s;
end
%%%REMOVAL of The segment the problem areaa
for k = 0:N-nw
sub = 0;
for n = 0:N-1
sub = sub + x(n+1).*w(n+1).*exp(-2i * pi * n * k / N);
end
figure;
%T = 0:N/fs;
%imagesc(T,freq(1:N/2),20*log10(abs(output(1:N/2)')))
%axis xy; axis tight; colormap(jet);
%ylim([0 1000]);
%xlabel('Time in sec');
%ylabel('Frequency (Hz)');
plot(freq(1:N/2),20*log10(abs(output(1:N/2)')))
grid on

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!