Apply a time window to my fft

I have an excel file (to big to attach) that I am reading a signal from it is in time and volts. I want to apply and time window to my fft to make it increment along the signal and give me the frequencies. Now with my current MATLAB set up I do NOT have signal toolbox and other such tools, therefore I need to construct a for-loop to do my iterations a window if you will. I need help applying the for-loop. I have attached my current code, I believe that if this is done right I will end up with and array of frequencies and then an array of velocities that I can plot.

I believe that the for-loop should go before the fft and end before the %plot frequency.

Answer 1

[Star Strider]

My version of MATLAB cannot run images of code.

Perhaps something like this would work:

t = linspace(0, 100, 1000);
s = sin(21*pi*5*t);
Ts = t(2)-t(1);
Fs = 1/Ts;
Fn = Fs/2;
seglen = 100;
Fv = linspace(0, 1, fix(seglen/2)+1)*Fs;
Iv = 1:numel(Fv);
FTs = zeros(fix(numel(t)/seglen), seglen);
for k = 1:10
    idx = k : k+seglen-1;
    FTs(k,:) = fft(s(idx))/seglen;
    seg(k) = k;
end

figure
surf(Fv, seg, abs(FTs(:,Iv))*2)
grid on
xlabel('Frequency')
ylabel('Time (segment)')
title('STFT Matrix')

figure
surf(Fv, seg, abs(FTs(:,Iv))*2, 'EdgeColor','none')
grid on
xlabel('Frequency')
ylabel('Time (segment)')
view(90,90)
title('Spectrogram')
axis('tight')

Experiment to get appropriate results with your signals.

.

[shawn finley]

thanks for the help i am working on applying it. I did find a way to put my code in took some doing cuz i was using a virtual machine.

clc

clear all

close all

%variables

% the below data range ('A22:C250147') should be the Laser wavelenght

%(prior to event) should have no velocity i.e. v=0

signal=readtable('U296.xlsx','Range','A22:B250147'); %signal prior to event

class(signal); %returns 'table'

head(signal); % tells the variable names of the data sheet

laswaveL=((1528+1563)/2)*10^-9; % this is the average of the range given units m

Fs=2.5e9; % sampling frequency in Hz

L=length(signal.x_0_003);

T=1/Fs;

tvec=(0:L-1)*T;

%plot raw data

figure(1)

time=signal.x_0_0002;

volts=signal.x_0_003;

plot(time,volts)

title('Raw signal')

xlabel('time (sec)')

ylabel('volts (v)')

grid on

%going form time domain to frequency domain wiht fft this is to get

%frequency that is used to calc velocity

n=2^nextpow2(L/2);

f=Fs*(0:(L/2))/L;

x=fft(signal.x_0_003,n)/L;

P2=abs(x/L);

P1=P2(1:L/2+1);

P1(2:end-1)=2*P1(2:end-1);

%plot frequencey

figure(2)

plot(f,P1)

xlabel('f (Hz)')

ylabel('PSD (Power Spectrum Density |P1(f)|)')

[~,loc]=max(P1);

freq=f(loc)% units Hz

title(['Frequency estimate',num2str(freq),'Hz'])

grid on

%velocity at current time interval

velocity=(laswaveL/(4*pi))*freq %unis m/s

%plot velocity

%figure(3)

%plot(tvec,velocity)

%xlabel('time (sec)')

%ylabel('velocity (m/s)')

%grid on

[Star Strider]

My pleasure!

I can’t run that, so I’ll assume it’s correct and produces the desired result. I don’t see anything that looks like my code it that.

Note —

This divides ‘x’ by ‘L^2’, not ‘L’:

x=fft(signal.x_0_003,n)/L;
P2=abs(x/L);

That could produce amplitudes lower that you might otherwise expect.

Also, the frequency vector calculation needs to consider ‘n’, not ‘L’ as the vector length, and it must go from 0 Hz (d-c) to the Nyquist frequency (‘Fn’), with ‘n/2’ elements for a one-sided Fourier transform, and ‘-Fn’ to ‘Fn’ with ‘n’ elements for a two-sided Fourier transform.

If my Answer helped you solve your problem, please Accept it!

.

[shawn finley]

I have adjusted my code but i am getting "Index exceeds array bounds" error I believe that this is due to not having the for k=##:## configured correctly.

signal=readtable('U296.xlsx','Range','A22:B250147'); %signal prior to event

class(signal); %returns 'table'

head(signal); % tells the variable names of the data sheet

laswaveL=((1528+1563)/2)*10^-9; % this is the average of the range given units m

Fs=2.5e9; % sampling frequency in Hz

L=length(signal.x_0_003); % signal length

Fn=Fs/2; %Nyquist frequency

Fv=linspace(0,1,fix(L/2)+1)*Fs; %frequency vector

Iv=1:numel(Fv); %indexing vector

%plot raw data

figure

time=signal.x_0_0002;

volts=signal.x_0_003;

plot(time,volts)

title('Raw signal')

xlabel('time (sec)')

ylabel('volts (v)')

grid on

%going form time domain to frequency domain wiht fft this is to get

%frequency that is used to calc velocity

n=2^nextpow2(L/2);

f=Fs*(0:(L/2))/L;

for k=1:L+1;

idx=k:k+L-1;

FTs(k,:)=fft(signal.x_0_003(idx))/L;

seg(k)=k;

end

figure

surf(Fv,seg,abs(FTs(:,Iv))*2)

grid on

xlabel('Frequency')

ylabel('Time (segment)')

title('STFT Matrix')

[Star Strider]

I obviously can’t run that.

What line is throwing the error?

The only thing that is immediately obvious is:

% for k=1:L+1 

should probably be:

% for k=1:L

instead, although my code uses the desired number of segments (10) not the vector length to determine the index.

Since I have no idea what you’re doing, I’ll just let it go at that.

A slight improvement on my code would be:

t = linspace(0, 100, 1000).';                           % Transpose To Column Vector
s = sin(21*pi*5*t);
Ts = t(2)-t(1);
Fs = 1/Ts;
Fn = Fs/2;
seglen = 100;
nseg = fix(numel(t)/seglen);                            % Added
Fv = linspace(0, 1, fix(seglen/2)+1)*Fs;
Iv = 1:numel(Fv);
FTs = zeros(seglen, nseg);
for k = 1:nseg                                          % Changed
    idx = k : k+seglen-1;
    FTs(:,k) = fft(s(idx))/seglen;
    seg(k) = k;
end

figure
surf(seg, Fv, abs(FTs(Iv,:))*2)
grid on
xlabel('Time (segment)')
ylabel('Frequency')
title('STFT Matrix')

figure
surf(seg, Fv, abs(FTs(Iv,:))*2, 'EdgeColor','none')
grid on
xlabel('Time (segment)')
ylabel('Frequency')
view(0,90)
title('Spectrogram')
axis('tight')

This changes the orientation of ‘t’ to a column vector, and everything else necessarily changes to be compatible with it. The code is actually a bit simpler and easier to understand as the result of the changes.

.

[shawn finley]

line 108 is throing the error (see below). I have tried the k=1:L but it returns the same error same line. I have shrunk my data and attached it along with my code.

clc

clear all

close all

%variables

% the below data range ('A22:C250147') should be the Laser wavelenght

%(prior to event) should have no velocity i.e. v=0

signal=readtable('U296_10000.xlsx','Range','A22:B10022');

class(signal); %returns 'table'

head(signal); % tells the variable names of the data sheet

laswaveL=((1528+1563)/2)*10^-9; % this is the average of the range given units m

Fs=2.5e9; % sampling frequency in Hz

L=length(signal.x_0_003); % signal length

Fn=Fs/2; %Nyquist frequency

Fv=linspace(0,1,fix(L/2)+1)*Fs; %frequency vector

Iv=1:numel(Fv); %indexing vector

%plot raw data

figure

time=signal.x_0_0002;

volts=signal.x_0_003;

plot(time,volts)

title('Raw signal')

xlabel('time (sec)')

ylabel('volts (v)')

grid on

%going form time domain to frequency domain wiht fft this is to get

%frequency that is used to calc velocity

n=2^nextpow2(L/2);

f=Fs*(0:(L/2))/L;

for k=1:L;

idx=k:k+L-1;

FTs(k,:)=fft(signal.x_0_003(idx))/L;

seg(k)=k;

end

figure

surf(Fv,seg,abs(FTs(:,Iv))*2)

grid on

xlabel('Frequency')

ylabel('Time (segment)')

title('STFT Matrix')

[Star Strider]

Your signal and my code gave this result —

uz = unzip('https://www.mathworks.com/matlabcentral/answers/uploaded_files/714127/U296_10000.zip')

uz = 1×1 cell array
    {'U296_10000.xlsx'}

T1 = readtable(uz{1}, 'VariableNamingRule','preserve')

T1 = 10001×3 table
       TIME        CH4 (PDV)     CH1 (Det)
    ___________    __________    _________

        -0.0002        -0.003         0.03
        -0.0002    -0.0026016    0.0099219
        -0.0002       -0.0026         0.01
        -0.0002       -0.0026         0.03
        -0.0002    -0.0030016    0.0099219
        -0.0002    -0.0034016    0.0099219
        -0.0002       -0.0034         0.03
        -0.0002    -0.0026016    -0.010078
        -0.0002       -0.0022        -0.01
        -0.0002    -0.0026016    0.0099219
        -0.0002    -0.0022016    -0.010078
        -0.0002       -0.0026        -0.01
        -0.0002    -0.0030016    -0.010078
    -0.00019999    -0.0026016    -0.010078
    -0.00019999        -0.003         0.01
    -0.00019999       -0.0026        -0.01

% t = linspace(0, 100, 1000).';                           % Transpose To Column Vector
% s = sin(21*pi*5*t);
t = T1.TIME;
s = T1.('CH4 (PDV)');
Ts = mean(diff(t))

Ts = 4.0000e-10

% Tsd = std(diff(t))
Fs = 1/Ts

Fs = 2.5000e+09

% Ts = t(2)-t(1);
% Fs = 1/Ts;
Fn = Fs/2;
seglen = 100;
nseg = fix(numel(t)/seglen);
Fv = linspace(0, 1, fix(seglen/2)+1)*Fs;
Iv = 1:numel(Fv);
FTs = zeros(seglen, nseg);
for k = 1:nseg
    idx = k : k+seglen-1;
    sv = s(idx) - mean(s(idx));                         % Added
%     FTs(:,k) = fft(s(idx))/seglen;
    FTs(:,k) = fft(sv)/seglen;                          % Added
    seg(k) = k;
end

figure
surf(seg, Fv, abs(FTs(Iv,:))*2)
colormap(turbo)
grid on
xlabel('Time (segment)')
ylabel('Frequency')
title('STFT Matrix')

figure
surf(seg, Fv, abs(FTs(Iv,:))*2, 'EdgeColor','none')
colormap(turbo)
grid on
xlabel('Time (segment)')
ylabel('Frequency')
view(0,90)
title('Spectrogram')
axis('tight')

The one modification was to subtract the mean of each segment from the segment. An alternate approach would be to subtract the mean of the entire vector from the entire vector, then remove the corresponding lines in the loop and replace them with tthe original line (currently commented-out).

The Signal Processing Toolbox pspectrum function would be preferable here.

Make appropriate changes to get the result you want.

.

[shawn finley]

Thanks!! very much i will impliment it in my code.

Why are you using segleg=100 and not the lengthe of the signal?

[shawn finley]

I would love to use the pspectrum but i do not have the Signal Processing Toolbox

[shawn finley]

what is seglen and why is it set to 100?

[Star Strider]

My pleasure!

The ‘seglen’ variable is simply the segment length for the time-windowed fft. It can be anything you want, providing that it is less than the signal length (preferably no more than half for a 2-segment fft). The length of the fft itself can also be whatever you want, so even for a segment length of 100, the fft length can be 128, 256, 1024 or 4096 or any power of 2 (ideally).

Certain other changes would then be necessary, for example:

uz = unzip('https://www.mathworks.com/matlabcentral/answers/uploaded_files/714127/U296_10000.zip')

uz = 1×1 cell array
    {'U296_10000.xlsx'}

T1 = readtable(uz{1}, 'VariableNamingRule','preserve')

T1 = 10001×3 table
       TIME        CH4 (PDV)     CH1 (Det)
    ___________    __________    _________

        -0.0002        -0.003         0.03
        -0.0002    -0.0026016    0.0099219
        -0.0002       -0.0026         0.01
        -0.0002       -0.0026         0.03
        -0.0002    -0.0030016    0.0099219
        -0.0002    -0.0034016    0.0099219
        -0.0002       -0.0034         0.03
        -0.0002    -0.0026016    -0.010078
        -0.0002       -0.0022        -0.01
        -0.0002    -0.0026016    0.0099219
        -0.0002    -0.0022016    -0.010078
        -0.0002       -0.0026        -0.01
        -0.0002    -0.0030016    -0.010078
    -0.00019999    -0.0026016    -0.010078
    -0.00019999        -0.003         0.01
    -0.00019999       -0.0026        -0.01

% t = linspace(0, 100, 1000).';                           % Transpose To Column Vector
% s = sin(21*pi*5*t);
t = T1.TIME;
s = T1.('CH4 (PDV)');
Ts = mean(diff(t))

Ts = 4.0000e-10

% Tsd = std(diff(t))
Fs = 1/Ts

Fs = 2.5000e+09

% Ts = t(2)-t(1);
% Fs = 1/Ts;
Fn = Fs/2;
seglen = 75;                                            % Shorten 'seglen' to 75 (Changed)
nseg = fix(numel(t)/seglen);
nfft = 4096;                                            % FFT Length (Added)
Fv = linspace(0, 1, fix(nfft/2)+1)*Fs;
Iv = 1:numel(Fv);
FTs = zeros(nfft, nseg);
for k = 1:nseg
    idx = k : k+seglen-1;
    sv = s(idx) - mean(s(idx));                         % Added
%     FTs(:,k) = fft(s(idx))/seglen;
    FTs(:,k) = fft(sv,nfft)/seglen;                     % Added - Changed
    seg(k) = k;
end

figure
surf(seg, Fv, abs(FTs(Iv,:))*2)
colormap(turbo)
grid on
xlabel('Time (segment)')
ylabel('Frequency')
title('STFT Matrix')

figure
surf(seg, Fv, abs(FTs(Iv,:))*2, 'EdgeColor','none')
colormap(turbo)
grid on
xlabel('Time (segment)')
ylabel('Frequency')
view(0,90)
title('Spectrogram')
axis('tight')

This results in somewhat smoother fft results. With the latest changes here, ‘nfft’ can be anything you want (within the memory capacity of the machine it’s being run on), and the rest of the code will adapt.

.

[shawn finley]

thanks i will keep that in mind as I make the code do what I want (or at least what I think I want).

One last quetion I have been seeing when I look for examples of this type of analysis the below code could you try and explain it

L=length(X) %which is length of signal

nfft=2^nextpow2(L/2) %which I assuem is the 128,256,ect.

x=fft(X,nfft)/L

f=Fs/2*linspace(0,1,nfft/2+1)

[Star Strider]

That simply produces the normalised Fourier transform.

The first line is obvious.

The second calculates the length of the fft as the power-of-2 greater than ‘L’. (That makes the fft calculation faster and more efficient, and slightly increases the frequency resolution.)

The third computes the Fourier transform and normalises it by length. (Those results would have to be multiplied by 2 to approximate the amplitude of the original waveform at those frequencies.)

The fourth calculates the associated frequency vector.

The only convenient vector that’s missing creates an index vector that is the length of ‘f’. It’s not required, however it’s easier than using ‘1:numel(f)’ everywhere to make ‘x’ and ‘f’ the same lengths, for plotting and certain other uses.

.

[shawn finley]

How dose seglen relate back to the time in the signal?

[Star Strider]

The ‘siglen’ value does not relate to time. It is simply the vector signal length for the Fourier transform computation, and can be anything (within reason).

The time would actually be represented by the ‘seg’ vector. Here , it is the index in the loop, however it could be converted to time as:

seg(k) = median(t(idx));

or the first or last elements in the vector, the mean value of the time vector, or something else.

So in the context of the code:

uz = unzip('https://www.mathworks.com/matlabcentral/answers/uploaded_files/714127/U296_10000.zip')

uz = 1×1 cell array
    {'U296_10000.xlsx'}

T1 = readtable(uz{1}, 'VariableNamingRule','preserve')

T1 = 10001×3 table
       TIME        CH4 (PDV)     CH1 (Det)
    ___________    __________    _________

        -0.0002        -0.003         0.03
        -0.0002    -0.0026016    0.0099219
        -0.0002       -0.0026         0.01
        -0.0002       -0.0026         0.03
        -0.0002    -0.0030016    0.0099219
        -0.0002    -0.0034016    0.0099219
        -0.0002       -0.0034         0.03
        -0.0002    -0.0026016    -0.010078
        -0.0002       -0.0022        -0.01
        -0.0002    -0.0026016    0.0099219
        -0.0002    -0.0022016    -0.010078
        -0.0002       -0.0026        -0.01
        -0.0002    -0.0030016    -0.010078
    -0.00019999    -0.0026016    -0.010078
    -0.00019999        -0.003         0.01
    -0.00019999       -0.0026        -0.01

% t = linspace(0, 10, 250).';                           % Transpose To Column Vector
% s = sin(2*pi*5*t);
t = T1.TIME;
s = T1.('CH4 (PDV)');
Ts = mean(diff(t))

Ts = 4.0000e-10

% Tsd = std(diff(t))
Fs = 1/Ts

Fs = 2.5000e+09

% Ts = t(2)-t(1);
% Fs = 1/Ts;
Fn = Fs/2;
seglen = 50;                                            % Shorten 'seglen' to 50 (Changed)
nseg = fix(numel(t)/seglen);
nfft = 1024;                                            % FFT Length (Added)
Fv = linspace(0, 1, fix(nfft/2)+1)*Fn;
Iv = 1:numel(Fv);
FTs = zeros(nfft, nseg);
for k = 1:nseg
    idx = k : k+seglen-1;
    sv = s(idx) - mean(s(idx));                         % Added
%     FTs(:,k) = fft(s(idx))/seglen;
    FTs(:,k) = fft(sv,nfft)/seglen;                     % Added - Changed
    seg(k) = median(t(idx));
end

figure
hsp = surfc(seg, Fv, abs(FTs(Iv,:))*2);
colormap(turbo)
hsp(1).EdgeColor = 'none';
grid on
xlabel('Time (segment)')
ylabel('Frequency')
title('STFT Matrix')

figure
surf(seg, Fv, abs(FTs(Iv,:))*2, 'EdgeColor','none')
colormap(turbo)
grid on
xlabel('Time (segment)')
ylabel('Frequency')
view(0,90)
title('Spectrogram')
axis('tight')
Ax = gca;
Ax.TickDir = 'both';
Ax.XMinorTick = 'on';
Ax.YMinorTick = 'on';
Ax.XTickLabelRotation = 45;

figure
hcf = contourf(seg, Fv, abs(FTs(Iv,:))*2);
colormap(turbo)
grid on
xlabel('Time (segment)')
ylabel('Frequency')
title(' Contour Plot Of STFT Matrix')
Ax = gca;
Ax.TickDir = 'both';
Ax.XMinorTick = 'on';
Ax.YMinorTick = 'on';
Ax.XTickLabelRotation = 45;

I improved it a bit, and just for fun, added the contourf call.

Tweak it further to get different results.

.

[shawn finley]

Thank you fot the clairifacation.

If i wanted to calculate the PSD of the fft I would want to do that inside the loop correct?

[Star Strider]

As always, my pleasure.

‘If i wanted to calculate the PSD of the fft I would want to do that inside the loop correct?’

Yes. This code just calculates the Fourier transform, so one approach to calculating the PSD would be to simply square it (using element-wise operations, so .^), or express it as:

PSD = (abs(FTs(Iv,:))*2).^2

or:

PSD = 20*log10(abs(FTs(Iv,:))*2)

That would work with my existing code.

A more typical approach is the pwelch calculation, so one option is adapting my code to work with it:

uz = unzip('https://www.mathworks.com/matlabcentral/answers/uploaded_files/714127/U296_10000.zip')

uz = 1×1 cell array
    {'U296_10000.xlsx'}

T1 = readtable(uz{1}, 'VariableNamingRule','preserve')

T1 = 10001×3 table
       TIME        CH4 (PDV)     CH1 (Det)
    ___________    __________    _________

        -0.0002        -0.003         0.03
        -0.0002    -0.0026016    0.0099219
        -0.0002       -0.0026         0.01
        -0.0002       -0.0026         0.03
        -0.0002    -0.0030016    0.0099219
        -0.0002    -0.0034016    0.0099219
        -0.0002       -0.0034         0.03
        -0.0002    -0.0026016    -0.010078
        -0.0002       -0.0022        -0.01
        -0.0002    -0.0026016    0.0099219
        -0.0002    -0.0022016    -0.010078
        -0.0002       -0.0026        -0.01
        -0.0002    -0.0030016    -0.010078
    -0.00019999    -0.0026016    -0.010078
    -0.00019999        -0.003         0.01
    -0.00019999       -0.0026        -0.01

% t = linspace(0, 10, 250).';                           % Transpose To Column Vector
% s = sin(2*pi*5*t);
t = T1.TIME;
s = T1.('CH4 (PDV)');
Ts = mean(diff(t))

Ts = 4.0000e-10

% Tsd = std(diff(t))
Fs = 1/Ts

Fs = 2.5000e+09

% Ts = t(2)-t(1);
% Fs = 1/Ts;
Fn = Fs/2;
seglen = 50;                                            % Shorten 'seglen' to 50 (Changed)
nseg = fix(numel(t)/seglen);
nfft = 1024;                                            % FFT Length (Added)
Fv = linspace(0, 1, fix(nfft/2)+1)*Fn;
Iv = 1:numel(Fv);
% FTs = zeros(nfft, nseg);
FTs = zeros(2^nextpow2(seglen)*2+1, nseg);              % Change To Use 'pwelch'
f = linspace(0, 1, size(FTs,1))*Fn;                     % Apparently, 'pwelch' Cannot Work With 'Fs' As Defined Here
for k = 1:nseg
    idx = k : k+seglen-1;
    sv = s(idx) - mean(s(idx));                         % Added
%     FTs(:,k) = fft(s(idx))/seglen;
%     FTs(:,k) = fft(sv,nfft)/seglen;                     % Added - Changed
    FTs(:,k) = pwelch(sv);                              % Call 'pwelch'
    seg(k) = median(t(idx));
end

figure
% hsp = surfc(seg, Fv, abs(FTs(Iv,:))*2);
hsp = surfc(seg, f, FTs);
colormap(turbo)
hsp(1).EdgeColor = 'none';
grid on
xlabel('Time (segment)')
ylabel('Frequency')
title('STFT Matrix')

figure
surf(seg, f, FTs, 'EdgeColor','none')
colormap(turbo)
grid on
xlabel('Time (segment)')
ylabel('Frequency')
view(0,90)
title('Spectrogram')
axis('tight')
Ax = gca;
Ax.TickDir = 'both';
Ax.XMinorTick = 'on';
Ax.YMinorTick = 'on';
Ax.XTickLabelRotation = 45;

figure
hcf = contourf(seg, f, FTs);
colormap(turbo)
grid on
xlabel('Time (segment)')
ylabel('Frequency')
title(' Contour Plot Of STFT Matrix')
Ax = gca;
Ax.TickDir = 'both';
Ax.XMinorTick = 'on';
Ax.YMinorTick = 'on';
Ax.XTickLabelRotation = 45;

This simply adapts pwelch to my existing code, although changing any of the existing parameters may break it.. A more robust implementation is in the spectrogram and similar functions.

.