Merge branch 'develop'

chanvocoder.m

@@ -0,0 +1,57 @@
+function y = chanvocoder(carrier, modul, chan, numband, overlap)
+% y = chanvocoder(carrier, modul, chan, numband, overlap)
+% The Channel Vocoder modulates the carrier signal with the modulation signal
+% chan = number of channels         (e.g., 512)
+% numband = number of bands (<chan) (e.g., 32)
+% overlap = window overlap          (e.g., 1/4)
+
+if numband>chan
+	error('# bands must be < # channels')
+end
+
+[rc, cc]   = size(carrier);
+if cc>rc
+	carrier = carrier';
+end
+
+[rm, cm]   = size(modul);
+if cm>rm
+	modul = modul';
+end
+
+st         = min(rc,cc);                         % stereo or mono?
+if st~= min(rm,cm)
+	error('carrier and modulator must have same number of tracks');
+end
+
+len        = min(length(carrier),length(modul)); % find shortest length
+carrier    = carrier(1:len,1:st);                % shorten carrier if needed
+modul      = modul(1:len,1:st);                  % shorten modulator if needed
+L          = 2*chan;                             % window length/FFT length
+w          = hanning(L); 
+if st==2
+	w=[w w];
+end  % window/ stereo window
+
+bands      = 1:round(chan/numband):chan;         % indices for frequency bands     
+bands(end) = chan;
+y          = zeros(len,st);                      % output vector
+
+ii         = 0;
+while ii*L*overlap+L <= len
+    ind    = round([1+ii*L*overlap:ii*L*overlap+L]);
+    FFTmod = fft( modul(ind,:) .* w );    % window & take FFT of modulator
+    FFTcar = fft( carrier(ind,:) .* w );  % window & take FFT of carrier  
+    syn    = zeros(chan,st);              % place for synthesized output
+    for jj = 1:numband-1                  % for each frequency band
+        b        = [bands(jj):bands(jj+1)-1]; % current band
+        syn(b,:) = FFTcar(b,:)*diag(mean(abs(FFTmod(b,:))));
+    end                                   % take product of spectra
+    midval   = FFTmod(1+L/2,:).*FFTcar(1+L/2,:); % midpoint is special
+    synfull  = [syn; midval; flipud( conj( syn(2:end,:) ) );]; % + and - frequencies
+    timsig   = real( ifft(synfull) );     % invert back to time
+    y(ind,:) = y(ind,:) + timsig;         % add back into time waveform   
+    ii       = ii+1;
+end
+y = 0.8*y/max(max(abs(y)));               % normalize output
+

frequencySpectrum.m

@@ -0,0 +1,66 @@
+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)')
+
+hold off
+figure;
+plot(f,10*log10(power/power(ind)));
+xlabel('Frequency (Hz)')
+ylabel('Power (dB)')
+

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)');

speech_analysis.m

@@ -0,0 +1,217 @@
+function speech_analysis()
+clear all
+close all
+clc
+
+
+% Construct the full file path
+filepath = './sound/modulator22.wav';
+% Read the audio file
+[y, Fs] = audioread(filepath);
+disp(['Successfully read the audio file: ', filepath]);
+disp(['Sampling frequency (Fs): ', num2str(Fs), ' Hz']);
+disp(['Number of samples: ', num2str(length(y))]);
+
+% Construct the output filename correctly
+%[~, name, ~] = fileparts(filepath); % Get the filename without extension
+%outputFilename = fullfile('.', ['processed_', name, '.wav']); % Create the new filename
+
+% Write the audio to a new file with double the sampling rate
+%audiowrite(outputFilename, y, Fs*2);
+%disp(['Successfully wrote the processed audio to: ', outputFilename, ' with double the sampling rate.']);
+
+disp('Playing the audio with double the sampling rate.');
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Plot
+t = (0:length(y)-1) / Fs; % Time in seconds
+
+figure;
+plot(t, y);
+xlabel('Time (seconds)');
+ylabel('Amplitude');
+title(['Temporal Variation of ', filepath]);
+grid on;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Frequency Spectrum
+%FFT
+tic;
+[yFFT, FFT_Time]=frequencySpectrum(y,Fs, 1);
+disp(FFT_Time);
+%DFT
+tic
+[yDFT, DFT_Time]=frequencySpectrum(y,Fs, 0);
+disp(DFT_Time);
+
+%Modify the padding to make the change.
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% Spectrogram
+spectrogram(y, Fs, 5,50)
+title('Spectrogram of modulator22.wav');
+colorbar;
+ylabel('Frequency (Hz)');
+xlabel('Time (s)');
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%Going to Pratt, we see that :
+%F0 : (100 + 130 +  100 + 120 + 100 + 90) / 6
+%F1 : 578.3725189859462, 418.70239431349677, 552.8090680139439,
+%308.88658136343446, 314.17710770594937, 363.8180262223959
+%F2 : 1695.8136433413672, 1550.9109531347972, 566.7831612330604,
+%1721.8044733141373, 1802.7920754749957, 1891.9059418088873
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% First downsampling (Shannon-Nyquist problem)
+desiredFreq = 4000; %in Hz
+
+% --- Downsampling using downsample() ---
+downsample_factor_ds = round(Fs / desiredFreq);
+y_downsampled_ds = downsample(y, downsample_factor_ds);
+Fs_downsampled_ds = Fs / downsample_factor_ds;
+disp(['--- Downsampling using downsample() ---']);
+disp(['New sampling frequency (downsample): ', num2str(Fs_downsampled_ds), ' Hz']);
+disp(['Number of samples (downsample): ', num2str(length(y_downsampled_ds))]);
+
+% --- Downsampling using decimate() ---
+downsample_factor_dec = round(Fs / desiredFreq);
+y_decimated = decimate(y, downsample_factor_dec);
+Fs_decimated = Fs / downsample_factor_dec;
+disp(['--- Downsampling using decimate() ---']);
+disp(['New sampling frequency (decimate): ', num2str(Fs_decimated), ' Hz']);
+disp(['Number of samples (decimate): ', num2str(length(y_decimated))]);
+
+%% --- Plotting Downsampled Signals ---
+figure;
+subplot(3,1,1);
+plot(t, y);
+xlabel('Time (seconds)');
+ylabel('Amplitude');
+title(['Original Signal (Fs = ', num2str(Fs), ' Hz)']);
+grid on;
+t_ds = (0:length(y_downsampled_ds)-1) / Fs_downsampled_ds;
+subplot(3,1,2);
+plot(t_ds, y_downsampled_ds);
+xlabel('Time (seconds)');
+ylabel('Amplitude');
+title(['Downsampled Signal (downsample, Fs = ', num2str(Fs_downsampled_ds), ' Hz)']);
+grid on;
+t_dec = (0:length(y_decimated)-1) / Fs_decimated;
+subplot(3,1,3);
+plot(t_dec, y_decimated);
+xlabel('Time (seconds)');
+ylabel('Amplitude');
+title(['Decimated Signal (decimate, Fs = ', num2str(Fs_decimated), ' Hz)']);
+grid on;
+%{
+%% --- Frequency Spectrum of Downsampled Signals ---
+figure;
+subplot(2,1,1);
+[yFFT_ds, FFT_Time_ds]=frequencySpectrum(y_downsampled_ds,Fs_downsampled_ds, 1);
+disp(['FFT Time (downsampled): ', num2str(FFT_Time_ds)]);
+plot(yFFT_ds, Fs_downsampled_ds);
+title('FFT of Downsampled Signal (downsample)');
+subplot(2,1,2);
+[yFFT_dec, FFT_Time_dec]=frequencySpectrum(y_decimated,Fs_decimated, 1);
+disp(['FFT Time (decimated): ', num2str(FFT_Time_dec)]);
+plot(yFFT_dec, Fs_decimated)
+title('FFT of Decimated Signal (decimate)');
+%}
+%{
+% --- Spectrograms of Downsampled Signals ---
+figure;
+subplot(2,1,1);
+spectrogram(y_downsampled_ds, round(0.02*Fs_downsampled_ds), round(0.01*Fs_downsampled_ds), 512, Fs_downsampled_ds, 'yaxis');
+title(['Spectrogram of Downsampled Signal (downsample, Fs = ', num2str(Fs_downsampled_ds), ' Hz)']);
+colorbar;
+ylabel('Frequency (Hz)');
+xlabel('Time (s)');
+subplot(2,1,2);
+spectrogram(y_decimated, round(0.02*Fs_decimated), round(0.01*Fs_decimated), 512, Fs_decimated, 'yaxis');
+title(['Spectrogram of Decimated Signal (decimate, Fs = ', num2str(Fs_decimated), ' Hz)']);
+colorbar;
+ylabel('Frequency (Hz)');
+xlabel('Time (s)');
+%}
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+desiredFilterFreq = 1000;
+%% --- Low-pass FIR filter ---
+order_fir = 30;
+normalized_cutoff_fir = desiredFilterFreq / (Fs / 2);
+y_fir_coeffs = fir1(order_fir, normalized_cutoff_fir, 'low'); % 'low' specifies a low-pass filter
+y_fir_filtered = filter(y_fir_coeffs, 1, y); % Apply the FIR filter
+
+figure;
+freqz(y_fir_coeffs, 1, 512, Fs); % Plot the frequency response of the FIR filter
+title('Frequency Response of FIR Low-Pass Filter');
+
+% FIR filter stability check (always stable)
+disp('--- FIR Filter Stability ---');
+disp('FIR filters designed using fir1 are inherently stable.');
+
+% --- Low-pass IIR filter (Butterworth) ---
+order_iir = 8;
+normalized_cutoff_iir = desiredFilterFreq / (Fs / 2);
+[b_iir, a_iir] = butter(order_iir, normalized_cutoff_iir, 'low'); % 'low' specifies a low-pass filter
+y_iir_filtered = filter(b_iir, a_iir, y); % Apply the IIR filter
+
+figure;
+freqz(b_iir, a_iir, 512, Fs); % Plot the frequency response of the IIR filter
+title('Frequency Response of IIR (Butterworth) Low-Pass Filter');
+%{
+% IIR filter stability check
+disp('--- IIR Filter (Butterworth) Stability ---');
+poles_iir = roots(a_iir);
+magnitudes_iir = abs(poles_iir);
+if all(magnitudes_iir < 1)
+    disp('The IIR (Butterworth) filter is stable (all poles are inside the unit circle).');
+else
+    disp('The IIR (Butterworth) filter is NOT stable (some poles are outside or on the unit circle).');
+    disp('Poles magnitudes:');
+    disp(magnitudes_iir);
+end
+%}
+
+%% --- Downsampling after filtering ---
+downsample_factor_filtered = round(Fs / desiredFreq);
+Fs_ds_filtered = Fs / downsample_factor_filtered;
+%{
+y_ds_fir_filtered = downsample(y_fir_filtered, downsample_factor_filtered);
+disp(['--- Downsampling FIR filtered signal using downsample() ---']);
+disp(['New sampling frequency (FIR filtered, downsample): ', num2str(Fs_ds_filtered), ' Hz']);
+disp(['Number of samples (FIR filtered, downsample): ', num2str(length(y_ds_fir_filtered))]);
+%}
+y_ds_iir_filtered = downsample(y_iir_filtered, downsample_factor_filtered);
+disp(['--- Downsampling IIR filtered signal using downsample() ---']);
+disp(['New sampling frequency (IIR filtered, downsample): ', num2str(Fs_ds_filtered), ' Hz']);
+disp(['Number of samples (IIR filtered, downsample): ', num2str(length(y_ds_iir_filtered))]);
+
+%% Plotting the signals
+% --- Comparing Output Signals ---
+% Temporal Variation
+figure;
+subplot(3,1,1);
+plot(t, y);
+xlabel('Time (seconds)');
+ylabel('Amplitude');
+title('Original Signal');
+grid on;
+subplot(3,1,2);
+plot(t, y_fir_filtered);
+xlabel('Time (seconds)');
+ylabel('Amplitude');
+title('FIR Filtered Signal');
+grid on;
+subplot(3,1,3);
+plot(t, y_iir_filtered);
+xlabel('Time (seconds)');
+ylabel('Amplitude');
+title('IIR Filtered Signal');
+grid on;
+
+% Play audios (using the audio data 'y' and its sampling rate 'Fs')
+%sound(y, Fs); % Play the original sound
+%sound(y, Fs*2);
+%sound(y_decimated,Fs_decimated)
+%sound(y_downsampled_ds,Fs_downsampled_ds) %Has distortion. This is because the Shannon-Nyquist criteria is not respected. Downsample() doesn't make sure the signal is filtered. Decimate does. So if need to choose, choose decimate !
+
+end