From 80e5f9a983d6aef6ad4508f001bfaf0c0fa49aa4 Mon Sep 17 00:00:00 2001 From: Loic Delattre Date: Fri, 24 Mar 2023 08:13:58 +0100 Subject: [PATCH] initial function files --- chanvocoder.m | 57 ++++++++++++++++++++++++++++++++++++++++++ frequencySpectrum.m | 60 +++++++++++++++++++++++++++++++++++++++++++++ spectrogram.m | 38 ++++++++++++++++++++++++++++ 3 files changed, 155 insertions(+) create mode 100644 chanvocoder.m create mode 100644 frequencySpectrum.m create mode 100644 spectrogram.m diff --git a/chanvocoder.m b/chanvocoder.m new file mode 100644 index 0000000..09fad17 --- /dev/null +++ b/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 + 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 + diff --git a/frequencySpectrum.m b/frequencySpectrum.m new file mode 100644 index 0000000..6d3a67a --- /dev/null +++ b/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)') + diff --git a/spectrogram.m b/spectrogram.m new file mode 100644 index 0000000..f01e8b5 --- /dev/null +++ b/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