-
Notifications
You must be signed in to change notification settings - Fork 0
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
2 changed files
with
139 additions
and
139 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,62 +1,62 @@ | ||
| function [x, k] = istft(S, window, hop, symFlag) | ||
| % Inverse Short Time Fourier Transform in the least squares sense. | ||
| % Given an STFT array S, this function finds x such that | ||
| % norm( S - stft(x, window, hop, size(S,1)) , 'fro' ) | ||
| % is minimized. | ||
| % | ||
| % The algorithm is based on [1]. | ||
| % The implementation is fully vectorized and supports multichannel signals. | ||
| % | ||
| % | ||
| % Inputs: | ||
| % S - STFT array, with frequency along rows and time across columns. | ||
| % If size(S,3) > 1, the ISTFT is performed independently on each layer. | ||
| % window - a vector containing the window function. | ||
| % hop - a positive integer describing the hop size. | ||
| % symFlag (optional) - frequency symmetry type, specified as 'nonsymmetric' (default) or 'symmetric'. Same as symFlag of MATLAB's built-in ifft. | ||
| % Set as 'symmetric' if x is expected to be real. | ||
| % Outputs: | ||
| % x - a time domain signal. For 3D S, the the i'th column of x if the ISTFT of S(:,:,i). | ||
| % k - a scalar describing the instability of the inverse transform. It depends only on the window and hop size. | ||
| % The minimum value is 1, which means the most stable. | ||
| % k=1 is achived for windows and hop sizes whos square value satisfy the COLA (Constant Overlap Add) condition. | ||
| % | ||
| % [1] Griffin, Daniel, and Jae Lim. "Signal estimation from modified short-time Fourier transform." IEEE Transactions on Acoustics, Speech, and Signal Processing 32.2 (1984): 236-243. | ||
| % | ||
| % Written by Tom Shlomo, 2019 | ||
| %% defauls | ||
| if nargin<4 | ||
| symFlag = 'nonsymmetric'; | ||
| end | ||
|
|
||
| %% Name some dimensions | ||
| N = size(S,2); | ||
| L = length(window); | ||
| Q = size(S,3); | ||
| T = L + hop*(N-1); | ||
|
|
||
| %% IFFT to obtain segemnts | ||
| y = ifft(S, [], 1, symFlag); | ||
| y(L+1:end,:,:) = []; % this is for case where NFFT is larger than the window's length | ||
|
|
||
| %% Weighted overlap add | ||
| y = y.*window(:); % apply weights | ||
| I = (1:L)' + hop*(0:N-1) + reshape( (0:Q-1)*T, 1, 1, [] ); % calculate subs for accumarray | ||
| x = reshape( accumarray(I(:), y(:)), [], Q ); % overlap-add | ||
|
|
||
| %% Least-squares correction | ||
| L2 = ceil(L/hop)*hop; | ||
| window(end+1:L2) = 0; % zero pad window to be an integer multiple of hop | ||
| J = mod( (1:hop:L2)' + (0:hop-1) -1, L2 ) + 1; | ||
| denom = sum(window(J).^2, 1)'; | ||
| if nargout==2 | ||
| k = max(denom)/min(denom); | ||
| end | ||
| denom = repmat(denom, ceil(T/hop), 1); | ||
| denom(T+1:end) = []; | ||
| x = x./denom; | ||
|
|
||
| %% delete the first L-1 rows of x | ||
| x(1:L-1, :) = []; | ||
|
|
||
| function [x, k] = istft(S, window, hop, symFlag) | ||
| % Inverse Short Time Fourier Transform in the least squares sense. | ||
| % Given an STFT array S, this function finds x such that | ||
| % norm( S - stft(x, window, hop, size(S,1)) , 'fro' ) | ||
| % is minimized. | ||
| % | ||
| % The algorithm is based on [1]. | ||
| % The implementation is fully vectorized and supports multichannel signals. | ||
| % | ||
| % | ||
| % Inputs: | ||
| % S - STFT array, with frequency along rows and time across columns. | ||
| % If size(S,3) > 1, the ISTFT is performed independently on each layer. | ||
| % window - a vector containing the window function. | ||
| % hop - a positive integer describing the hop size. | ||
| % symFlag (optional) - frequency symmetry type, specified as 'nonsymmetric' (default) or 'symmetric'. Same as symFlag of MATLAB's built-in ifft. | ||
| % Set as 'symmetric' if x is expected to be real. | ||
| % Outputs: | ||
| % x - a time domain signal. For 3D S, the the i'th column of x if the ISTFT of S(:,:,i). | ||
| % k - a scalar describing the instability of the inverse transform. It depends only on the window and hop size. | ||
| % The minimum value is 1, which means the most stable. | ||
| % k=1 is achived for windows whos square value satisfy the COLA (Constant Overlap Add) condition. | ||
| % | ||
| % [1] Griffin, Daniel, and Jae Lim. "Signal estimation from modified short-time Fourier transform." IEEE Transactions on Acoustics, Speech, and Signal Processing 32.2 (1984): 236-243. | ||
| % | ||
| % Written by Tom Shlomo, 2019 | ||
| %% defauls | ||
| if nargin<4 | ||
| symFlag = 'nonsymmetric'; | ||
| end | ||
|
|
||
| %% Name some dimensions | ||
| N = size(S,2); | ||
| L = length(window); | ||
| Q = size(S,3); | ||
| T = L + hop*(N-1); | ||
|
|
||
| %% IFFT to obtain segemnts | ||
| y = ifft(S, [], 1, symFlag); | ||
| y(L+1:end,:,:) = []; % this is for case where NFFT is larger than the window's length | ||
|
|
||
| %% Weighted overlap add | ||
| y = y.*window(:); % apply weights | ||
| I = (1:L)' + hop*(0:N-1) + reshape( (0:Q-1)*T, 1, 1, [] ); % calculate subs for accumarray | ||
| x = reshape( accumarray(I(:), y(:)), [], Q ); % overlap-add | ||
|
|
||
| %% Least-squares correction | ||
| L2 = ceil(L/hop)*hop; | ||
| window(end+1:L2) = 0; % zero pad window to be an integer multiple of hop | ||
| J = mod( (1:hop:L2)' + (0:hop-1) -1, L2 ) + 1; | ||
| denom = sum(window(J).^2, 1)'; | ||
| if nargout==2 | ||
| k = max(denom)/min(denom); | ||
| end | ||
| denom = repmat(denom, ceil(T/hop), 1); | ||
| denom(T+1:end) = []; | ||
| x = x./denom; | ||
|
|
||
| %% delete the first L-1 rows of x | ||
| x(1:L-1, :) = []; | ||
|
|
||
| end |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -1,79 +1,79 @@ | ||
| function [S, f, t] = stft(x, window, hop, NFFT, fs) | ||
| % Short Time Fourier Transform of the signal x. | ||
| % The implementation is fully vectorized and supports multichannel signals. | ||
| % | ||
| % Inputs: | ||
| % x The time domain signal. must be either a column vector or a 2D array. | ||
| % If x is a 2D array, the STFT is calculated for each column independently. | ||
| % window A vector containing the window function. | ||
| % Default value is hann(128,'periodic'). | ||
| % hop hop size | ||
| % NFFT Number of DFT points. | ||
| % Default value is 2^nextpow2(length(window)) | ||
| % fs Sample frequency. Doesn't affect S, only f and t. | ||
| % Default value is 1. | ||
| % Outputs: | ||
| % S STFT array, with frequency across rows, and time across columns. | ||
| % If x has multiple columms, than S(:,:,i) is the STFT of x(:,i). | ||
| % f Frequency vector. Same unit as fs. | ||
| % t Time vector (starting from 0). Same units as 1/fs. | ||
| % | ||
| % Please note that unlike MATLAB's builtin function spectrogram, this | ||
| % implementation of the STFT is causal. For example, if hop=1, then S(:,i) corresponds to a | ||
| % segment of x whos *last* index is i. More precisely: | ||
| % S(:,i,j) = fft(x(i-L+1:i,j).*window, NFFT) | ||
| % Where L = length(window), and x is zeropadded where neccessary. | ||
| % This defnition has several advantages over the more standard, anti-causal one. | ||
| % First, the inverse transform is more stable near the begining of the signal. | ||
| % Second, each frequency of the STFT can be regarded as the (resampled) | ||
| % output of an FRI filter. | ||
| % | ||
| % | ||
| % Written by Tom Shlomo, 2019 | ||
|
|
||
|
|
||
| %% default values | ||
| if nargin<2 || isempty(window) | ||
| window = hann(128,'periodic'); | ||
| end | ||
| windowLen = length(window); | ||
|
|
||
| if nargin<3 || isempty(hop) | ||
| hop = max(1, round(windowLen/4)); | ||
| end | ||
|
|
||
| if nargin<4 || isempty(NFFT) | ||
| NFFT = 2^nextpow2(windowLen); | ||
| end | ||
|
|
||
| if nargin<5 || isempty(fs) | ||
| fs = 1; | ||
| end | ||
|
|
||
| %% Padd with zeros from the left (for causality) | ||
| x = [zeros(windowLen-1,size(x,2)); x]; | ||
|
|
||
| %% partition x to data segments | ||
| n = size(x,1); | ||
| I = 1:hop:n; | ||
| m = (0:windowLen-1)'; | ||
| I = I+m; | ||
| x(end+1:I(end), :)=0; | ||
| I = I + reshape( (0:size(x,2)-1)*size(x,1) , 1, 1, [] ); | ||
| x = x(I); | ||
|
|
||
| %% apply window function | ||
| x = x.*window(:); | ||
|
|
||
| %% FFT each segment | ||
| S = fft(x, NFFT, 1); | ||
|
|
||
| %% calculate frequency and time vectors | ||
| if nargout>=2 | ||
| f = (0:NFFT-1)'/NFFT * fs; | ||
| end | ||
| if nargout>=3 | ||
| t = (0:size(S,2)-1)'*hop/fs; | ||
| end | ||
|
|
||
| function [S, f, t] = stft(x, window, hop, NFFT, fs) | ||
| % Short Time Fourier Transform of the signal x. | ||
| % The implementation is fully vectorized and supports multichannel signals. | ||
| % | ||
| % Inputs: | ||
| % x The time domain signal. must be either a column vector or a 2D array. | ||
| % If x is a 2D array, the STFT is calculated for each column independently. | ||
| % window A vector containing the window function. | ||
| % Default value is hann(128,'periodic'). | ||
| % hop hop size | ||
| % NFFT Number of DFT points. | ||
| % Default value is 2^nextpow2(length(window)) | ||
| % fs Sample frequency. Doesn't affect S, only f and t. | ||
| % Default value is 1. | ||
| % Outputs: | ||
| % S STFT array, with frequency across rows, and time across columns. | ||
| % If x has multiple columms, than S(:,:,i) is the STFT of x(:,i). | ||
| % f Frequency vector. Same unit as fs. | ||
| % t Time vector (starting from 0). Same units as 1/fs. | ||
| % | ||
| % Please note that unlike MATLAB's builtin function spectrogram, this | ||
| % implementation of the STFT is causal. For example, if hop=1, then S(:,i) corresponds to a | ||
| % segment of x whos *last* index is i. More precisely: | ||
| % S(:,i,j) = fft(x(i-L+1:i,j).*window, NFFT) | ||
| % Where L = length(window), and x is zeropadded where neccessary. | ||
| % This defnition has several advantages over the more standard, anti-causal one. | ||
| % First, the inverse transform is more stable near the begining of the signal. | ||
| % Second, each frequency of the STFT can be regarded as the (resampled) | ||
| % output of an FIR filter. | ||
| % | ||
| % | ||
| % Written by Tom Shlomo, 2019 | ||
|
|
||
|
|
||
| %% default values | ||
| if nargin<2 || isempty(window) | ||
| window = hann(128,'periodic'); | ||
| end | ||
| windowLen = length(window); | ||
|
|
||
| if nargin<3 || isempty(hop) | ||
| hop = max(1, round(windowLen/4)); | ||
| end | ||
|
|
||
| if nargin<4 || isempty(NFFT) | ||
| NFFT = 2^nextpow2(windowLen); | ||
| end | ||
|
|
||
| if nargin<5 || isempty(fs) | ||
| fs = 1; | ||
| end | ||
|
|
||
| %% Padd with zeros from the left (for causality) | ||
| x = [zeros(windowLen-1,size(x,2)); x]; | ||
|
|
||
| %% partition x to data segments | ||
| n = size(x,1); | ||
| I = 1:hop:n; | ||
| m = (0:windowLen-1)'; | ||
| I = I+m; | ||
| x(end+1:I(end), :)=0; | ||
| I = I + reshape( (0:size(x,2)-1)*size(x,1) , 1, 1, [] ); | ||
| x = x(I); | ||
|
|
||
| %% apply window function | ||
| x = x.*window(:); | ||
|
|
||
| %% FFT each segment | ||
| S = fft(x, NFFT, 1); | ||
|
|
||
| %% calculate frequency and time vectors | ||
| if nargout>=2 | ||
| f = (0:NFFT-1)'/NFFT * fs; | ||
| end | ||
| if nargout>=3 | ||
| t = (0:size(S,2)-1)'*hop/fs; | ||
| end | ||
|
|
||
| end |