MATLAB: Is it possible to find the frequency specture from a signal imported from data points (excel)

disturbancefrequency

I measured a signal with a lot of disturbance. The data was exported to excel and i had imported the data in Matlab and made a plot. Now i want to create a frequency plot. (where i can see the power of each frequency, to see which frequencies are wrong in the signal. I tried to use fft and ifft but with no results.

Do you have suggestions?

Thanks in advance

Lucas

Best Answer

hello

for your help see my little demo below for audio wav file analysis. Just adapt the load to your data file format

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [signal,Fs]=wavread('audio_in_noise1.wav');  %(older matlab)
% or 
[data,Fs]=audioread('audio_in_noise1.wav'); %(newer matlab)
channel = 1;
signal = data(:,channel);
samples = length(signal);
dt = 1/Fs;
t = (0:dt:(samples-1)*dt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
samples = length(signal);
NFFT = 512*4;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum, freq] = pwelch(signal,w,NOVERLAP,NFFT,Fs);
figure(1),semilogx(freq,20*log10(sensor_spectrum));grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
saturation_dB = 80;  % dB range scale (means , the lowest displayed level is 80 dB below the max level)
[sg,fsg,tsg] = specgram(signal,NFFT,Fs,w,NOVERLAP);  
% 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
% saturation to given dB range
mini_dB = round(max(max(sg_dBpeak))) - saturation_dB;
sg_dBpeak(sg_dBpeak<mini_dB) = mini_dB;
% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);axis('xy');colorbar('vert');grid
title(['Spectrogram  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');

Related Solutions

MATLAB: Speech Signal Analysis using IIR filters

hello

your code is about frequency analysis (FFT) and not filtering time series. ... so what is the purpose of bandpass filtering here ?

what info are you looking at ?

maybe a time frequency analysis can be of some help to understand energy distribution vs time and vs frequency. Then maybe after that first analysis you can filter the signal for another use.

If it can help you, below some code to do basic spectral and time / frequency analysis. I can further help you on the digital filters once we understand what you are looking for.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% dummy signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fs = 5000;
samples = 255001;
dt = 1/Fs;
t = (0:dt:(samples-1)*dt);
omega = 2*pi*(25+20*t);
sensordata = randn(size(t)) + 5*sin(omega.*t);  % signal = linear chirp + noise
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
samples = length(sensordata);
NFFT = 512;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% option 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum, freq] = pwelch(sensordata,w,NOVERLAP,NFFT,Fs);
figure(1),plot(freq,20*log10(sensor_spectrum));
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Amplitude (dB)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% option 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
saturation_dB = 80;  % dB range scale (means , the lowest displayed level is 80 dB below the max level)
% fmin = 1;
% fmax = Fs/2;
[sg,fsg,tsg] = specgram(sensordata,NFFT,Fs,w,NOVERLAP);  
% 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
% saturation to given dB range
mini_dB = round(max(max(sg_dBpeak))) - saturation_dB;
sg_dBpeak(sg_dBpeak<mini_dB) = mini_dB;
% plots spectrogram
figure(2);
imagesc(tsg,fsg,sg_dBpeak);axis('xy');colorbar('vert');
title(['Spectrogram  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');

MATLAB: How to only pass the signals that I want

hello

this code will do fft / spectrogram and bandpass filtering - adapt the corner frequencies to your own needs !

clc
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 4096;    % 
OVERLAP = 0.5;
% 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 = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% data
% [data,Fs] = audioread('myWAVaudiofile.wav');
[data,Fs] = audioread('Immigran_Hiss.wav');
channel = 1;
signal = data(:,channel);
samples = length(signal);
%% decimate (if needed)
% NB : decim = 1 will do nothing (output = input)
decim = 1;
if decim>1
    signal = decimate(signal,decim);
    Fs = Fs/decim;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% display 1 : 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(1),plot(freq,sensor_spectrum_dB,'b');grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel(my_ylabel);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



% display 2 : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sg,fsg,tsg] = specgram(signal,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);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% removal of hiss between t = 0 and t = 1.68 s
% use of bandpass filter to remove noise below 500 and above 2500 Hz
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data_filtered = data;
t_min = 0;
t_max = 1.68;
ind = 1+fix(t_min*Fs:t_max*Fs-1);
% keep only signal content extracted by bandpass filters 
f_low = 500;
f_high = 2500;
N = 4;
% signal_filtered = zeros(size(signal));
[b,a] = butter(N,2/Fs*[f_low f_high]);
% now let's filter the signal
data_filtered(ind,:) = filtfilt(b,a,data(ind,:));
% save as wav file 
audiowrite('Immigran_Hiss-filt.wav',data_filtered,Fs);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 2B : time / frequency analysis : spectrogram demo
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sg,fsg,tsg] = specgram(data_filtered(:,1),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(3);
imagesc(tsg,fsg,sg_dBpeak);colormap('jet');
axis('xy');colorbar('vert');grid
title([my_title ' / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(fsg(2)-fsg(1)) ' Hz ']);
xlabel('Time (s)');ylabel('Frequency (Hz)');
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 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;
    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

Best Answer

Related Solutions

MATLAB: Speech Signal Analysis using IIR filters

MATLAB: How to only pass the signals that I want

Related Question