extend signal_from_spectrum to nd matrices

Author: Fiete Winter

SFS_general/signal_from_spectrum.m

@@ -1,4 +1,4 @@
-function outsig = signal_from_spectrum(amplitude,phase,f,conf)
+function outsig = signal_from_spectrum(amplitude,phase,f,dim,conf)
 %SIGNAL_FROM_SPECTRUM time signal from single-sided spectrum
 %
 %   Usage: outsig = signal_from_spectrum(amplitude,phase,f,conf)
@@ -7,6 +7,7 @@
 %       amplitude   - the single-sided amplitude spectrum
 %       phase       - the single-sided phase spectrum / rad
 %       f           - the corresponding frequency vector
+%       dim         - dimension along which the fft is performed
 %       conf        - configuration struct (see SFS_config)
 %
 %   Output parameters:
@@ -50,10 +51,24 @@
 
 %% ===== Checking input arguments ========================================
 nargmin = 4;
-nargmax = 4;
+nargmax = 5;
 narginchk(nargmin,nargmax);
-[amplitude,phase] = column_vector(amplitude,phase);
+if nargin == nargmin
+   conf = dim;
+   dim = find(size(amplitude)~=1, 1);  % find first non-singleton dimension 
+else   
+   isargpositivescalar(dim); 
+end
+isargstruct(conf);
 
+Ndims = ndims(amplitude);
+dim = min(dim,Ndims);
+amplitude = permute(amplitude, [dim:Ndims, 1:dim-1]);  % move dim to first dimension
+phase = permute(phase, [dim:Ndims, 1:dim-1]);  % move dim to first dimension
+s = size(amplitude);
+Nx = s(1);
+amplitude = reshape(amplitude, Nx, []);  % squeeze all other dimensions
+phase = reshape(phase, Nx, []);  % squeeze all other dimensions
 
 %% ===== Configuration ===================================================
 fs = conf.fs;
@@ -67,25 +82,31 @@
     % Length of the signal to generate
     samples = 2 * (bins-1);
     % Rescaling (see spectrum_from_signal())
-    amplitude = [amplitude(1); amplitude(2:end-1)/2; amplitude(end)] * samples;
+    amplitude = [amplitude(1,:); amplitude(2:end-1,:)/2; amplitude(end,:)] ...
+        * samples;
     % Mirror the amplitude spectrum ( 2*pi periodic [0, fs[ )
-    amplitude = [amplitude; amplitude(end-1:-1:2)];
+    amplitude = [amplitude; amplitude(end-1:-1:2,:)];
     % Mirror the phase spectrum and build the inverse (complex conjugate)
-    phase = [phase; -1 * phase(end-1:-1:2)];
+    phase = [phase; -1 * phase(end-1:-1:2,:)];
 
 else  % -> odd time signal length
     % Length of the signal to generate
     samples = 2*bins - 1;
     % Rescaling (see signal_from_spectrum)
-    amplitude = [amplitude(1); amplitude(2:end)/2] * samples;
+    amplitude = [amplitude(1,:); amplitude(2:end,:)/2] * samples;
     % Mirror the amplitude spectrum ( 2*pi periodic [0, fs-bin] )
-    amplitude = [amplitude; amplitude(end:-1:2)];
+    amplitude = [amplitude; amplitude(end:-1:2,:)];
     % Mirror the phase spectrum and build the inverse (complex conjugate)
-    phase = [phase; -1*phase(end:-1:2)];
+    phase = [phase; -1*phase(end:-1:2,:)];
 end
 
 % Convert to complex spectrum
 compspec = amplitude .* exp(1i*phase);
 
 % Build the inverse fft and assume spectrum is conjugate symmetric
 outsig = real(ifft(compspec));
+
+%% ===== Output ==========================================================
+% undo reshape and permute
+outsig = reshape(outsig, [samples, s(2:end)]);
+outsig = permute(outsig, [Ndims-dim+2:Ndims, 1:Ndims-dim+1]); 

SFS_general/spectrum_from_signal.m

@@ -1,10 +1,11 @@
-function varargout = spectrum_from_signal(signal,conf)
+function varargout = spectrum_from_signal(signal,dim,conf)
 %SPECTRUM_FROM_SIGNAL single-sided amplitude and phase spectra of signal
 %
-%   Usage: [amplitude,phase,f] = spectrum_from_signal(sig,conf)
+%   Usage: [amplitude,phase,f] = spectrum_from_signal(sig,[dim],conf)
 %
 %   Input parameters:
 %       signal      - one channel audio (time) signal
+%       dim         - dimension along which the fft is performed
 %       conf        - configuration struct (see SFS_config)
 %
 %   Output parameters:
@@ -52,11 +53,22 @@
 
 %% ===== Check input arguments ===========================================
 nargmin = 2;
-nargmax = 2;
+nargmax = 3;
 narginchk(nargmin,nargmax);
-signal = column_vector(signal);
+if nargin == nargmin
+   conf = dim;
+   dim = find(size(signal)~=1, 1);  % find first non-singleton dimension 
+else   
+   isargpositivescalar(dim); 
+end
 isargstruct(conf);
 
+Ndims = ndims(signal);
+dim = min(dim,Ndims);
+signal = permute(signal, [dim:Ndims, 1:dim-1]);  % move dim to first dimension
+s = size(signal);
+Nx = s(1);
+signal = reshape(signal, Nx, []);  % squeeze all other dimensions
 
 %% ===== Configuration ===================================================
 fs = conf.fs;
@@ -65,7 +77,7 @@
 
 %% ===== Calcualate spectrum =============================================
 % Generate fast fourier transformation (=> complex output)
-compspec = fft(signal);
+compspec = fft(signal,[],1);
 
 % Length of the signal => number of points of fft
 bins = length(signal);
@@ -75,37 +87,47 @@
     f = fs/bins * (0:(bins-1)/2)';
     % Get amplitude and phase spectra and use only the first half of the
     %>spectrum [0, fs/2[
-    amplitude = abs(compspec(1:length(f)));
-    phase = angle(compspec(1:length(f)));
+    amplitude = abs(compspec(1:length(f),:));
+    phase = angle(compspec(1:length(f),:));
     % Scale the amplitude (factor two for mirrored frequencies
     %>divide by number of bins)
-    amplitude = [amplitude(1); 2*amplitude(2:end)] / bins;
+    amplitude = [amplitude(1,:); 2*amplitude(2:end,:)] / bins;
 
 else  % For even signal length
     % Calculate corresponding frequency axis
     f = fs/bins * (0:bins / 2)';
     % Get amplitude and phase spectra and use only the first half of the
     %>spectrum [0, fs/2]
-    amplitude = abs(compspec(1:length(f)));
-    phase = angle(compspec(1:length(f)));
+    amplitude = abs(compspec(1:length(f),:));
+    phase = angle(compspec(1:length(f),:));
     % Scale the amplitude (factor two for mirrored frequencies
     %>divide by number of bins)
-    amplitude = [amplitude(1); 2*amplitude(2:end-1); amplitude(end)] / bins;
+    amplitude = [amplitude(1,:); 2*amplitude(2:end-1,:); amplitude(end,:)] / bins;
 end
 
-% Return values
-if nargout>0, varargout{1}=amplitude; end
-if nargout>1, varargout{2}=phase; end
-if nargout>2, varargout{3}=f; end
-
-
 %% ===== Plotting ========================================================
 if nargout==0 || useplot
     figure; title('Spectrum');
     subplot(2,1,1)
     semilogx(f,20 * log10(abs(amplitude))); xlim([1, fs/2]);
     grid on; xlabel('frequency / Hz'); ylabel('amplitude / dB')
     subplot(2,1,2)
-    semilogx(f,unwrap(phase)); xlim([1, fs/2]);
+    semilogx(f,unwrap(phase,[],1)); xlim([1, fs/2]);
     grid on; xlabel('frequency / Hz'); ylabel('phase / rad')
 end
+
+%% ===== Output ==========================================================
+% Return values
+if nargout>0
+    % undo reshape and permute
+    amplitude = reshape(amplitude, [length(f), s(2:end)]);
+    amplitude = permute(amplitude, [Ndims-dim+2:Ndims, 1:Ndims-dim+1]); 
+    varargout{1}=amplitude;
+end
+if nargout>1
+    % undo reshape and permute
+    phase = reshape(phase, [length(f), s(2:end)]);
+    phase = permute(phase, [Ndims-dim+2:Ndims, 1:Ndims-dim+1]);
+    varargout{2}=phase; 
+end
+if nargout>2, varargout{3}=f; end