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

5 次查看（过去 30 天）
Aniket Manjare
Aniket Manjare 2020-12-19
评论： Walter Roberson 2020-12-20
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

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

回答（1 个）

Mathieu NOE
Mathieu NOE 2020-12-20
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(:);

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by