Script for audio analyzing

12 views (last 30 days)
Incoronata Gagliardi
Incoronata Gagliardi on 9 Nov 2022
Commented: Mathieu NOE on 9 Nov 2022
Hello, i really hope someone can help me, basically here is my problem:
i'm working on a script for audio spectrum and peak search, we used an oscilloscope to export .csv data of a monocromatic audio wave for fft peak search. I want to read this file on matlab and plot the fft but i don't really know how to do it. I successfully imported the .csv file in which there are three columns, times, amplitude and fft trasform, but when i plot the 3rd one with times or amplitude it doesn't show me the peak. I tryed other ways to plot the signal like dividing the frequencies span and plotting the fft with frequences but nothing change and the peak doesn't appear. I'm really sorry if i missed some terms and wrote in a bad way, but still i need some help because i cannot find a way to plot the data.
Thank you if anyone will answer this.
  2 Comments
Jonas
Jonas on 9 Nov 2022
can you upload an exemplary csv, then we can have a look.
are fft values complex or real valued? you can start by just plotting the 3rd column, using plot(theColumn) or plot(abs(theColumn)) if values are complex
Incoronata Gagliardi
Incoronata Gagliardi on 9 Nov 2022
Here is the csv file, the peak should be a 80 Hz, i tried to plot only 3rd column but nothing changes, as you can see plotting it the result is a peak in 0 and nothing noticeable after that. Here is the base script i did.
Please forgive me if there are some mistakes, i know something is not working and i'm worried everything is wrong, i tried to "play" with codes, but nothing good comes from them. If you can help me, i will immensly grateful, but i'm sorry to disturb you so feel free to not help and anyway thank you :).
T=readtable('test peak 80 Hz0.csv');
% ^^^^^^^^^------ name csv file
t=T{:,1};
amp=T{:,2};
fft=T{:,3};
l_fft=length(fft)/2; %because only half of fft column are not NaN
frq_max=1000; %max frequency
frq_min=0; %min frequency
dfrq=(frq_max-frq_min)/(l_fft-1);
colonna4=frq_min:dfrq:frq_max; %here i divided the frequency span to have a set of "points"
colonna4=transpose(colonna4); %need the column
fft_estratta=fft([1:l_fft],:); %here i had some problems with dimensions of arrays, so i cut at the length i needed
colonna4=colonna4([1:l_fft],:); %needed only some values
plot(colonna4, fft_estratta, 'r')

Sign in to comment.

Answers (1)

Mathieu NOE
Mathieu NOE on 9 Nov 2022
hello
try this code
FYI, having the sampling rate at 1MHz is massive overkill if you intend to analyse only below 1 kHz. It makes just the data file 500 times too big.
Had to implement the decimation before fft to avoid a 1 million bins fft computation !
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
filename = ('test peak 80 Hz0.csv');
data = readmatrix(filename);
time = data(:,1);
signal = data(:,2);
dt = mean(diff(time));
Fs = 1/dt; % sampling rate
[samples,channels] = size(signal);
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 500;
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;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 1000; % so df = 2 Hz (Fs / NFFT)
OVERLAP = 0.75;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : time domain plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1),
plot(time,signal);grid on
title(['Time plot / Fs = ' num2str(Fs) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[freq, spectrum_raw] = myfft_peak(signal,Fs,NFFT,OVERLAP);
figure(2),plot(freq,20*log10(spectrum_raw));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('Amplitude (dB)');
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
  3 Comments
Incoronata Gagliardi
Incoronata Gagliardi on 9 Nov 2022
I really, really appreciate your help, i have no words to thank you. I'm going to test it tomorrow with my laboratory mates and will let you know if the problem is fixed, really thank you!!

Sign in to comment.

Products


Release

R2022b

Community Treasure Hunt

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

Start Hunting!