facing issue in FFT : P1 = P2(1:L/2+1)

14 views (last 30 days)
Aniket Manjare
Aniket Manjare on 19 Dec 2020
Commented: Walter Roberson on 20 Dec 2020
I am facing error to solve this code , kindly guide me
i am facing issue at:
P1 = P2(1:L/2+1)
Kindly help me resove the issue: the error details screenshot has be uploaded
Below is the full code i am trying
%%
%Data Collection and Labeling
clc;
clear all;
labelPoints = {[2,50],[51,80]};
NoOfFailureModes = 2;
estimateRULFlags = [false, true];
path = 'Datasets/';
labels = {[categorical("OFF"),categorical("ON"),categorical("DUCT BLOCKAGE")],[categorical("OFF"),categorical("ON"),categorical("ROTOR IMBALANCE")]};
rawdata = cell(NoOfFailureModes);
data = cell(NoOfFailureModes);
for i = 1:NoOfFailureModes
rawdata{i} = readtable([path 'Dataset' num2str(i) '.csv'],'Delimiter',','); %Read the downloaded csv file
data{i} = dataParsing(rawdata{i},labelPoints{i},labels{i},i); %Parse the rawdata
end
%%
%Data Exploration
fs = 100; %Sampling Frequency
L = 100; %Length of each data sample is 100 sensor readings
f = fs*(0:(L/2))/L; % Frequency bins until Nyquist Frequency at which FFT is computed
convertToMinutes = (1/fs)/60; %Convert sample number to minutes 0.01/60
votype = 'VideoWriter'; %Record the Exploration to view it again
vo = VideoWriter('Data Exploration Video', 'MPEG-4');
set(vo, 'FrameRate', 5);
open(vo);
for i = 1:NoOfFailureModes
NoOfDataPoints = height(data{i});
x = zeros(L*NoOfDataPoints,0);
y = zeros(L*NoOfDataPoints,0);
z = zeros(L*NoOfDataPoints,0);
figure;
set(gcf,'Visible','on')
for j = 1:NoOfDataPoints
data{i}.X_FFT(j,:) = computeFFT(data{i}.X(j,:)-mean(data{i}.X(j,:)),L);
data{i}.Y_FFT(j,:) = computeFFT(data{i}.Y(j,:)-mean(data{i}.Y(j,:)),L);
data{i}.Z_FFT(j,:) = computeFFT(data{i}.Z(j,:)-mean(data{i}.Z(j,:)),L);
x((j-1)*L+1 : j*L) = data{i}.X(j,:);
y((j-1)*L+1 : j*L) = data{i}.Y(j,:);
z((j-1)*L+1 : j*L) = data{i}.Z(j,:);
t = (1:j*L)*convertToMinutes;
subplot(3,2,1);
plot(t,x);
title({"Experiment: " + num2str(i) , "Time Domain"});
xlabel("time (min.)",'FontWeight','bold');
ylabel("X",'rotation',0,'FontWeight','bold');
subplot(3,2,2);
bar(f,data{i}.X_FFT(j,:));
title({"Sample: " + num2str(j) + "; Label: " + char(data{i}.Label(j)), "Frequency Domain"});
xlabel("frequency (Hz)",'FontWeight','bold');
ylabel("X",'rotation',0,'FontWeight','bold');
subplot(3,2,3);
plot(t,y);
xlabel("time (min.)",'FontWeight','bold');
ylabel("Y",'rotation',0,'FontWeight','bold');
subplot(3,2,4);
bar(f,data{i}.Y_FFT(j,:));
xlabel("frequency (Hz)",'FontWeight','bold');
ylabel("Y",'rotation',0,'FontWeight','bold');
subplot(3,2,5);
plot(t,z);
xlabel("time (min.)",'FontWeight','bold');
ylabel("Z",'rotation',0,'FontWeight','bold');
subplot(3,2,6);
bar(f,data{i}.Z_FFT(j,:));
xlabel("frequency (Hz)",'FontWeight','bold');
ylabel("Z",'rotation',0,'FontWeight','bold');
drawnow;
writeVideo(vo, getframe(gcf));
end
end
close(vo);
%%
%%
%Helper Functions
function P1 = computeFFT(X,L)
Y = fft(X);
P2 = abs(Y/L);%absolute value
P1 = P2(1:L/2+1)
P1(2:end-1) = 2*P1(2:end-1);
end
function data = dataParsing(rawdata,labelPoints,labels,dataset)
bias = 44;
noOfDataPoints = length(rawdata.entry_id);
labelPoints = [0 labelPoints noOfDataPoints];
for i = 1:length(labelPoints)-1
Label(labelPoints(i)+1:labelPoints(i+1),:) = labels(i); %Create Labels for the data points
end
j = 1;
if(dataset == 2)
offset = 0;
else
offset = 2000;
end
for i = 1:noOfDataPoints
var1 = rawdata.field1(i);
var2 = rawdata.field2(i);
var3 = rawdata.field3(i);
var4 = rawdata.field4(i);
var5 = rawdata.field5(i);
var6 = rawdata.field6(i);
%Turning the device off and on will reset the first data set to all zeroes
if var1 ~= 0 || var2 ~= 0 || var3 ~= 0
var = [var1, var2, var3, var4, var5, var6] - bias;
data.Label(j,:) = Label(i);
data.sno(j,:) = rawdata.entry_id(i) + offset;
data.X(j,:) = [var(1);var(4)];
data.Y(j,:) = [var(2);var(5)];
data.Z(j,:) = [var(3);var(6)];
j = j+1;
end
end
data = struct2table(data);
end

Answers (1)

Mathieu NOE
Mathieu NOE on 20 Dec 2020
Hi
suggest this fft code that also allows averaging if needed :
function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
%FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer overlap % (between 0 and 0.95)
samples = length(signal);
% fill signal with zeros if its length is lower than nfft
if samples<nfft
s_tmp = zeros(nfft,1);
s_tmp((1:samples)) = signal;
signal = s_tmp;
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;
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
freq_vector = freq_vector(:);
fft_spectrum = fft_spectrum(:);

Community Treasure Hunt

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

Start Hunting!

Translated by