117 lines
3.1 KiB
Matlab
117 lines
3.1 KiB
Matlab
function [power, duration] = frequencySpectrum(signal, fs, pad, verbose, rangefr)
|
|
%%%%%%%%%%%%%%%%%%
|
|
%function power = frequencySpectrum(signal, fs, pad)
|
|
%
|
|
% Task: Display the power spectrum (lin and log scale) of a given signal
|
|
%
|
|
% Input:
|
|
% - signal: the input signal to process
|
|
% - fs: the sampling rate
|
|
% -pad: boolean if true, signal is padded with 0 to the next power of 2 -> FFT instead of DFT
|
|
%
|
|
% Output:
|
|
% - power: the power spectrum
|
|
%
|
|
%
|
|
% Guillaume Gibert, guillaume.gibert@ecam.fr
|
|
% 25/04/2022
|
|
%%%%%%%%%%%%%%%%%%
|
|
|
|
n = length(signal); % number of samples
|
|
|
|
if (pad)
|
|
n = 2^nextpow2(n);
|
|
end
|
|
|
|
tic
|
|
y = fft(signal, n);% compute DFT of input signal
|
|
duration = toc;
|
|
|
|
power = abs(y).^2/n; % power of the DFT
|
|
|
|
[val, ind] = max(power); % find the mx value of DFT and its index
|
|
|
|
rangemin = rangefr(1);
|
|
rangemax = rangefr(2);
|
|
|
|
if (verbose)
|
|
% plots
|
|
figure;
|
|
|
|
subplot(1,3,1) % time plot
|
|
t=0:1/fs:(n-1)/fs; % time range
|
|
%pad signal with zeros
|
|
if (pad)
|
|
signal = [ signal; zeros( n-length(signal), 1)];
|
|
end
|
|
plot(t, signal)
|
|
xticks(0:0.1*fs:n*fs);
|
|
xticklabels(0:0.1:n/fs);
|
|
xlabel('Time (s)');
|
|
ylabel('Amplitude (a.u.)');
|
|
|
|
subplot(1,3,2) % linear frequency plot
|
|
f = (0:n-1)*(fs/n); % frequency range
|
|
plot(f,power, 'b*'); hold on;
|
|
plot(f,power, 'r');
|
|
xlabel('Frequency (Hz)')
|
|
ylabel('Power (a.u.)')
|
|
xlim([rangemin rangemax]);
|
|
|
|
|
|
subplot(1,3,3) % log frequency plot
|
|
power_db = 10*log10(power/power(ind));
|
|
plot(f, power_db);
|
|
xlabel('Frequency (Hz)')
|
|
ylabel('Power (dB)')
|
|
xlim([rangemin rangemax]);
|
|
|
|
subplot(1,3,3) % log frequency plot
|
|
power_db = 10*log10(power/power(ind));
|
|
plot(f, power_db);
|
|
xlabel('Frequency (Hz)')
|
|
ylabel('Power (dB)')
|
|
xlim([rangemin rangemax]);
|
|
|
|
% Shift power_db to all positive values
|
|
shift_amount = min(power_db) * -1;
|
|
power_db_shifted = power_db + shift_amount;
|
|
|
|
% Find peaks and their indices
|
|
[pks, locs] = findpeaks(power_db_shifted);
|
|
|
|
% Interpolate to get a smooth line
|
|
xi = linspace(min(f(locs)), max(f(locs)), 1000);
|
|
yi = interp1(f(locs), pks, xi, 'pchip');
|
|
|
|
% Shift the interpolated values back down and up by 10
|
|
yi = yi - shift_amount + 10;
|
|
|
|
% Smooth the line using a moving average filter
|
|
windowSize = 15; % Adjust this value to change the amount of smoothing
|
|
b = (1/windowSize)*ones(1,windowSize);
|
|
yi = filter(b, 1, yi);
|
|
yi = filter(b, 1, flip(yi)); % Apply the filter in reverse direction
|
|
yi = flip(yi); % Flip the data back to original direction
|
|
|
|
% Shift yi to all positive values
|
|
shift_amount_yi = min(yi) * -1;
|
|
yi_shifted = yi + shift_amount_yi;
|
|
|
|
% Find peaks of the smoothed line
|
|
[pks_smooth, locs_smooth] = findpeaks(yi_shifted);
|
|
|
|
% Plot the envelope
|
|
hold on;
|
|
plot(xi, yi, 'r-'); % plot the envelope
|
|
|
|
% Plot points and labels at the peaks
|
|
for i = 1:length(pks_smooth)
|
|
plot(xi(locs_smooth(i)), pks_smooth(i) - shift_amount_yi, 'ko');
|
|
text(xi(locs_smooth(i)), pks_smooth(i) - shift_amount_yi, ['F' num2str(i)], 'VerticalAlignment', 'bottom');
|
|
end
|
|
|
|
hold off;
|
|
|
|
end
|