FFT of large sample data

hey!

I have a very large data with 5M samples from my oscilloscope and im trying to make an fft of my data so it resemblance that of an spectrum analyzer with averaged value every 9k hz.

My approach was to calculate how many points i need to get a value every 9k hz over my spectrum from 0-50Mhz, whereas im actually interested in the spectrum between 150k to 30Mhz.

Then I made a matrix with the amount of points i wanted (11k) and then looped all the 5M samples into that matrix and took the average. When i measured a square wave with a known amplitude and frequency my code got the appropriate answer. BUT I wonder if the loop can actually be used like that when measuring real electrical circuits.

I would be thankful for any input!

 clear 
opts = detectImportOptions("flyback_nofilter_v1.csv");
opts.DataLines = [ 6 inf ];
V1 = readmatrix("flyback_nofilter_v1.csv",opts);


Ts = mean(diff(V1(:,1)));
Fs = 1/Ts;
Fn = Fs/2;
n = 1;
batchSize = 11111;
batches = length(V1) / batchSize;
StoredData = zeros(batchSize,1);
iterCount = 1;

 while iterCount < batches
     Data = V1(n:n+batchSize-1,2);
     StoredData = StoredData + Data;
     n = n + batchSize;
     iterCount = iterCount + 1;
 end

 Average = StoredData/batches;
LF = length(Average);
FT_v1 = ((fft(Average)*2)/LF); 
Fv = linspace(0,1,(LF/2)+1)*Fn;
Iv = 1:length(Fv);

dBuV_v1 = (20*log10(FT_v1)+120);

figure(1)
plot(Fv, (abs(dBuV_v1(Iv))));
line([150e3,500e3],[56,46])
line([500e3,5e6],[46,46])
line([5e6,30e6],[50,50])
line([5e6,5e6],[46,50])
 xlabel('Frequency (Hz)')
 ylabel('Amplitude (dBuV)')
set(gca, 'XLim',[150e3 30e6])
set(gca,'xscale','log')
 grid

[Mathieu NOE]

hello

there' something that disturbs me... you are averaging the time data , and then you do the fft on the averaged data ?

usually we do the averaging on the fft , unless you have a trigger signal that make you sure the time data buffer will overlay perfectly (to reduce the effect of spurious noise).

fyi , the averaging fft loop code :

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 percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    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*ones(1,channels));
        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

[Joakim Jennel]

Thanks for pointing that out, it is the main reason i made this post, to see if that is an acceptable solution or not. So yeah, right now Im averaging the amplitude in the time domain and afterwards Im doing the fft. Im getting the data samples from an Oscilloscope in a stop-state, im not using any trigger.

This was something that was suggested to me when I first started to write this code. At first i was thinking it sounded about right, you have X amount of curves that i just summed up together and averaged out and it would be fine. Later on when running the code on real measurements from different circuits I didnt get exactly the results i wanted. I was sure if it was because the code is incorrect or if it is simply to basic to get me accurate enough of a plot.

Debating whenever the average would be before or after the fft. Later on I was also thinking that it should be the otherway around but someone on this forum made a simple code showing that the mean value stayed the same, so I didnt really think that much about it. I was more leaning towards my loop in itself was incorrect.

In the first line of the code you sent "(signal, Fs, nfft, Overlap)" could the signal be a matrix of 5M data samples? And im not really sure what "overlap" is in my case. I have seen NFFT when I have been googling but Im not entirely sure how to use nfft in my case.

[Mathieu NOE]

hello

could you share some sample data ? I will show you how to use my code

all the best

[Joakim Jennel]

oh that would be great, thank you!

[Mathieu NOE]

ok , so here it is - avaraged spectrum (and spectrogram).

I was somehow assuming a kind of repetitive pattern but seems not to be the case. I have to admit it's not the usual data I'm processing (noise and vibration) so results may or not be meaningfull right now

this is my "work horse" code I'm using quite often :

clc
clearvars

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

filename = ('Data.csv'); 
data = readmatrix(filename,"NumHeaderLines",5); % 2 channels : Time, Ampl
time = data(:,1);
dt = mean(diff(time));
signal = data(:,2);
Fs = 1/dt; % sampling rate 

[samples,channels] = size(signal);

% time vector 
t = (0:samples-1)*dt;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
df = 9e3; % target delta f = 9kHz
NFFT = round(Fs/df);    % = 11111
% NFFT = 1024*4;    % 
OVERLAP = 0.75;

% spectrogram dB scale
spectrogram_dB_scale = 80;  % dB range scale (means , the lowest displayed level is XX dB below the max level)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% options 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% if you are dealing with acoustics, you may wish to have A weighted
% spectrums 

% option_w = 0 : linear spectrum (no weighting dB (L) )
% option_w = 1 : A weighted spectrum (dB (A) )
option_w = 0;

%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    for ck = 1:channels
    newsignal(:,ck) = decimate(signal(:,ck),decim);
    Fs = Fs/decim;
    end
   signal = newsignal;
end
samples = length(signal);
time = (0:samples-1)*1/Fs;

%%%%%% legend structure %%%%%%%%
for ck = 1:channels
    leg_str{ck} = ['Channel ' num2str(ck) ];
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),plot(time,signal);grid on
title(['Time plot  / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
legend(leg_str);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[freq, sensor_spectrum] = myfft_peak(signal,Fs,NFFT,OVERLAP);

% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum);

% apply A weigthing if needed
if option_w == 1
    pondA_dB = pondA_function(freq);
    sensor_spectrum_dB = sensor_spectrum_dB+pondA_dB;
    my_ylabel = ('Amplitude (dB (A))');
else
    my_ylabel = ('Amplitude (dB (L))');
end


figure(2),semilogx(freq,sensor_spectrum_dB);grid on
df = freq(2)-freq(1); % frequency resolution 
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
legend(leg_str);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 3 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for ck = 1:channels
    [sg,fsg,tsg] = specgram(signal(:,ck),NFFT,Fs,hanning(NFFT),floor(NFFT*OVERLAP));  

    % FFT normalisation and conversion amplitude from linear to dB (peak)
    sg_dBpeak = 20*log10(abs(sg))+20*log10(2/length(fsg));     % NB : X=fft(x.*hanning(N))*4/N; % hanning only

    % apply A weigthing if needed
    if option_w == 1
        pondA_dB = pondA_function(fsg);
        sg_dBpeak = sg_dBpeak+(pondA_dB*ones(1,size(sg_dBpeak,2)));
        my_title = ('Spectrogram (dB (A))');
    else
        my_title = ('Spectrogram (dB (L))');
    end

    % saturation of the dB range : 
    % saturation_dB = 60;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
    min_disp_dB = round(max(max(sg_dBpeak))) - spectrogram_dB_scale;
    sg_dBpeak(sg_dBpeak<min_disp_dB) = min_disp_dB;

    % plots spectrogram
    figure(2+ck);
    imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
    axis('xy');colorbar('vert');grid on
    df = fsg(2)-fsg(1); % freq resolution 
    title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df,3) ' Hz / Channel : ' num2str(ck)]);
    xlabel('Time (s)');ylabel('Frequency (Hz)');

end

function pondA_dB = pondA_function(f)
	% dB (A) weighting curve
	n = ((12200^2*f.^4)./((f.^2+20.6^2).*(f.^2+12200^2).*sqrt(f.^2+107.7^2).*sqrt(f.^2+737.9^2)));
	r = ((12200^2*1000.^4)./((1000.^2+20.6^2).*(1000.^2+12200^2).*sqrt(1000.^2+107.7^2).*sqrt(1000.^2+737.9^2))) * ones(size(f));
	pondA = n./r;
	pondA_dB = 20*log10(pondA(:));
end


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 percentage of overlap % (between 0 and 0.95)

[samples,channels] = size(signal);

% fill signal with zeros if its length is lower than nfft
if samples<nfft
    s_tmp = zeros(nfft,channels);
    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*ones(1,channels));
        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

[Joakim Jennel]

hey, thanks a lot! Now i have something to look at and in a better light.

Yeah, im measuring EMI noises and there are supposed to be spikes at various frequencies which I will later on try to fix with a filter. You usually measure this through a spectrum analyzer which I do not have, hence the simple code to get a better view at it.