how to eliminate zeros in the output of fast fourier transform (matrix / array)?

5 ビュー (過去 30 日間)
Naufal Arfani
Naufal Arfani 2021 年 1 月 4 日
コメント済み: Mathieu NOE 2021 年 1 月 5 日
I use simulink to find the frequency domain value of a sine signal that consists of two frequencies using fast fourier transform and then use findpeaks to find the value of each peak it gets. By using a zero order hold of 1e-3 and a buffer of 1000 and a time of 10 seconds so that the discrete matrix value becomes [11x1000] but there is an error so it can't be entered into findpeaks because I think it comes from the value t = 0 is all value 0 so that at [1x1000] there is a leading zero. I have used a lot of code but nothing works, is there any way to help me find out where the error is?
one example is I use the code below and I think because I omitted [1x1000] at the beginning so the rest becomes [10x1000] but it says an error because the number of elements doesn't match and I don't understand why.
in the way below is still the same, which part is lacking please I really need help
I also attach the workspace if possible it can be more helpful to understand
function y = fcn(u)
a = fftshift(u);
v = nonzeros(a');
y = reshape(v,1000,10);
function y = fcn(u)
y = u;
y(y==0) = [];
  3 件のコメント
Naufal Arfani
Naufal Arfani 2021 年 1 月 4 日
I'm sorry I mean the contents are like this so it looks like at t = 0 the value is all 0, but I tried to convert it to NaN but there were still errors

サインインしてコメントする。

採用された回答

Mathieu NOE
Mathieu NOE 2021 年 1 月 4 日
hello
FYI, this is a code to do averaged fft analysis + findpeaks demo + notch filter (extra bonus)
you can import the data from where you want
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% dummy data
Fs = 1000;
samples = 25000;
t = (0:samples-1)'*1/Fs;
signal = cos(2*pi*50*t)+cos(2*pi*100*t)+1e-3*rand(samples,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000; %
Overlap = 0.75;
w = hanning(NFFT); % Hanning window / Use the HANN function to get a Hanning window which has the first and last zero-weighted samples.
%% notch filter section %%%%%%
% H(s) = (s^2 + 1) / (s^2 + s/Q + 1)
fc = 50; % notch freq
wc = 2*pi*fc;
Q = 5; % adjust Q factor for wider (low Q) / narrower (high Q) notch
% at f = fc the filter has gain = 0
w0 = 2*pi*fc/Fs;
alpha = sin(w0)/(2*Q);
b0 = 1;
b1 = -2*cos(w0);
b2 = 1;
a0 = 1 + alpha;
a1 = -2*cos(w0);
a2 = 1 - alpha;
% analog notch (for info)
num1=[1/wc^2 0 1];
den1=[1/wc^2 1/(wc*Q) 1];
% digital notch (for info)
num1z=[b0 b1 b2];
den1z=[a0 a1 a2];
freq = linspace(fc-1,fc+1,200);
[g1,p1] = bode(num1,den1,2*pi*freq);
g1db = 20*log10(g1);
[gd1,pd1] = dbode(num1z,den1z,1/Fs,2*pi*freq);
gd1db = 20*log10(gd1);
figure(1);
plot(freq,g1db,'b',freq,gd1db,'+r');
title(' Notch: H(s) = (s^2 + 1) / (s^2 + s/Q + 1)');
legend('analog','digital');
xlabel('Fréquence (Hz)');
ylabel(' dB')
% now let's filter the signal
signal_filtered = filtfilt(num1z,den1z,signal);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display : averaged FFT spectrum before / after notch filter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq,fft_spectrum] = myfft_peak(signal, Fs, NFFT, Overlap);
sensor_spectrum_dB = 20*log10(fft_spectrum);% convert to dB scale (ref = 1)
% demo findpeaks
df = mean(diff(freq));
[pks,locs]= findpeaks(sensor_spectrum_dB,'SortStr','descend','NPeaks',2);
[freq,fft_spectrum_filtered] = myfft_peak(signal_filtered, Fs, NFFT, Overlap);
sensor_spectrum_filtered_dB = 20*log10(fft_spectrum_filtered);% convert to dB scale (ref = 1)
figure(2),semilogx(freq,sensor_spectrum_dB,'b',freq,sensor_spectrum_filtered_dB,'r');grid
title(['Averaged FFT Spectrum / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
legend('before notch filter','after notch filter');
xlabel('Frequency (Hz)');ylabel(' dB')
text(locs+.02,pks,num2str(freq(locs)))
function [freq_vector,fft_spectrum] = myfft_peak(signal, Fs, nfft, Overlap)
% FFT peak spectrum of signal (example sinus amplitude 1 = 0 dB after fft).
% Linear averaging
% signal - input signal,
% Fs - Sampling frequency (Hz).
% nfft - FFT window size
% Overlap - buffer overlap % (between 0 and 0.95)
signal = signal(:);
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;
samples = nfft;
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
  4 件のコメント
Naufal Arfani
Naufal Arfani 2021 年 1 月 5 日
sorry I ask for your help again @Mathieu NOE can you help make fft coding but with the input signal from simulink so maybe you don't need to use t and port because it will use t from simulink, can you help me?
Mathieu NOE
Mathieu NOE 2021 年 1 月 5 日
hello
I loaded the mat file but in don't understand what the data represents. I was waiting for a vector of same length as time vector. What are those 11 vectors that i plotted ? that has litte to do with a sine signal of 2 frequencies
>> load('file.mat')
>> out
out =
Simulink.SimulationOutput:
simout: [11x1000 double]
simout1: [1x1 double]
tout: [10001x1 double]
SimulationMetadata: [1x1 Simulink.SimulationMetadata]
ErrorMessage: [0x0 char]

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeSpectral Estimation についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by