Apply a time window to my fft
Show older comments
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.
Accepted Answer
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.
.
16 Comments
shawn finley
on 17 Aug 2021
Star Strider
on 17 Aug 2021
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
on 17 Aug 2021
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
on 17 Aug 2021
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')
Your signal and my code gave this result —
uz = unzip('https://www.mathworks.com/matlabcentral/answers/uploaded_files/714127/U296_10000.zip')
T1 = readtable(uz{1}, 'VariableNamingRule','preserve')
% 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))
% Tsd = std(diff(t))
Fs = 1/Ts
% 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).
Make appropriate changes to get the result you want.
.
shawn finley
on 17 Aug 2021
shawn finley
on 17 Aug 2021
shawn finley
on 17 Aug 2021
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')
T1 = readtable(uz{1}, 'VariableNamingRule','preserve')
% 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))
% Tsd = std(diff(t))
Fs = 1/Ts
% 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
on 17 Aug 2021
Star Strider
on 17 Aug 2021
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
on 18 Aug 2021
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')
T1 = readtable(uz{1}, 'VariableNamingRule','preserve')
% 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))
% Tsd = std(diff(t))
Fs = 1/Ts
% 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;
Tweak it further to get different results.
.
shawn finley
on 18 Aug 2021
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')
T1 = readtable(uz{1}, 'VariableNamingRule','preserve')
% 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))
% Tsd = std(diff(t))
Fs = 1/Ts
% 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.
.
More Answers (0)
Categories
Find more on Spectral Measurements in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
