Nice code

Signalanalysis.m

@@ -0,0 +1,44 @@
+clear all
+close all
+clc
+pkg load signal
+
+% Read the csv file
+signal = csvread('unknownsignal.csv');
+
+Fs = 300; % Sampling frequency
+t = (0:length(signal)-1)/Fs;
+N = length(signal);
+duration = N / Fs;
+% Plot the time-domain signal
+figure;
+plot(t, signal);
+xlabel('Time (s)');
+ylabel('Amplitude');
+title('Time-Domain Signal');
+window = hann(N)';
+windowed_data = window .* signal;
+a=size(windowed_data);
+display(a);
+% Perform spectral analysis using the spectrogram with specified parameters:
+
+
+[power,duration]=frequencySpectrum(signal, Fs, false);
+% Calculate the spectrogram
+
+
+
+
+% Design a Butterworth bandpass filter
+[b,a] = butter(3,0.4);
+y = filter(b, a, windowed_data);
+% Apply the filter to the signal
+spectrogram(y,Fs,5,30);
+
+% Visualize the filtered signal in the time domain:
+figure;
+plot(t, filtered_signal);
+xlabel('Time (s)');
+ylabel('Amplitude');
+title('Filtered Signal (30-40 Hz)');
+

frequencySpectrum.m

@@ -0,0 +1,60 @@
+function [power, duration] = frequencySpectrum(signal, fs, pad)
+%%%%%%%%%%%%%%%%%%
+%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
+
+% 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.)')
+
+subplot(1,3,3) % log frequency plot
+plot(f,10*log10(power/power(ind)));
+xlabel('Frequency (Hz)')
+ylabel('Power (dB)')
+

read.m

@@ -0,0 +1,42 @@
+clc
+clear all
+close all
+pkg load signal
+
+signal = csvread('unknownsignal.csv');
+Fs=300;
+N = length(signal);
+t = 0:1/Fs:(N-1)/Fs;
+rectwin = rectwin(N);
+figure;
+plot(t,signal.*rectwin);
+xlabel('Time (s)');
+ylabel('Amplitude');
+title('Unknown Signal in Time Domain');
+
+hammingwin = @hamming(N);
+signal_hamming = signal.*hammingwin;
+signal_fft = abs(fft(signal_hamming)/N).^2;
+f = Fs*(0:N-1)/N;
+figure;
+plot(f,signal_fft);
+xlabel('Frequency (Hz)');
+ylabel('Power Spectral Density (dB/Hz)');
+title('Power Spectral Density of Unknown Signal');
+
+Wn = [30 40]/(Fs/2);
+[b,a] = butter(4,Wn,'bandpass');
+signal_filtered = filter(b,a,signal);
+figure;
+plot(t,signal_filtered.*rectwin);
+xlabel('Time (s)');
+ylabel('Amplitude');
+title('Filtered Signal in Time Domain');
+
+signal_filtered_hamming = signal_filtered.*hammingwin;
+signal_filtered_fft = abs(fft(signal_filtered_hamming)/N).^2;
+figure;
+plot(f,signal_filtered_fft);
+xlabel('Frequency (Hz)');
+ylabel('Power Spectral Density (dB/Hz)');
+title('Power Spectral Density of Filtered Signal');

spectrogram.m

@@ -0,0 +1,38 @@
+function spectrogram(signal, samplingFreq, step_size, window_size)
+%%%%%%%%%%%%%%%%%%%%%%%
+%function spectrogram(signal, samplingFreq, step_size, window_size)
+% ex.:  spectrogram(signal, samplingFreq, step_size, window_size)
+%
+% Task: Plot the spectrogram of a given signal
+%
+% Inputs:
+%	-signal: temporal signal to analyse 
+%	-samplingFreq: sampling frequency of the temporal signal
+% 	-step_size: how often the power spectrum will be computed in ms
+%	-window_size: size of the analysing window in ms
+%
+% Ouput: None
+%
+% author: Guillaume Gibert (guillaume.gibert@ecam.fr)
+% date: 14/03/2023
+%%%%%%%%%%%%%%%%%%%%%%%
+
+figure;
+	subplot(2,1,1);
+t=0:1/samplingFreq:length(signal)/samplingFreq-1/samplingFreq;
+plot(t, signal');
+xlim([0 length(signal)/samplingFreq-1/samplingFreq]);
+ylabel('amplitude (norm. unit)');
+	subplot(2,1,2);
+step = fix(step_size*samplingFreq/1000);     % one spectral slice every step_size ms
+window = fix(window_size*samplingFreq/1000);  % window_size ms data window
+fftn = 2^nextpow2(window); % next highest power of 2
+[S, f, t] = specgram(signal, fftn, samplingFreq, window, window-step);
+S = abs(S(2:fftn*4000/samplingFreq,:)); % magnitude in range 0<f<=4000 Hz.
+S = S/max(S(:));           % normalize magnitude so that max is 0 dB.
+S = max(S, 10^(-40/10));   % clip below -40 dB.
+S = min(S, 10^(-3/10));    % clip above -3 dB.
+imagesc (t, f, log(S));    % display in log scale
+set (gca, "ydir", "normal"); % put the 'y' direction in the correct direction
+xlabel('time (s)');
+ylabel('frequency (Hz)');