Apply a time window to my fft

Question

shawn finley 2021-8-16

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/915627-apply-a-time-window-to-my-fft

评论： Star Strider 2021-8-18

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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2021-8-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/915627-apply-a-time-window-to-my-fft#answer_768352

在 MATLAB Online 中打开

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 个评论
显示 14更早的评论隐藏 14更早的评论

shawn finley 2021-8-17

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

shawn finley 2021-8-17

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 2021-8-17

在 MATLAB Online 中打开

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 2021-8-17

U296_10000.zip

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 2021-8-17

在 MATLAB Online 中打开

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.

.

Star Strider 2021-8-17

在 MATLAB Online 中打开

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.

.

Star Strider 2021-8-18

在 MATLAB Online 中打开

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 2021-8-18

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 2021-8-18

在 MATLAB Online 中打开

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.

.

请先登录，再进行评论。

Apply a time window to my fft

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

16 个评论
显示 14更早的评论隐藏 14更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Apply a time window to my fft

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

16 个评论 显示 14更早的评论隐藏 14更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

16 个评论
显示 14更早的评论隐藏 14更早的评论