Skip to content

Commit a9e1618

Browse files
dylewskydylewsky
authored
Add files via upload
1 parent 02dbef4 commit a9e1618

File tree

16 files changed

+428
-0
lines changed

16 files changed

+428
-0
lines changed

‎UTILS/CC2.mat‎

456 Bytes
Binary file not shown.

‎UTILS/SVT.m‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
function out = SVT(X,tau)
2+
[U,S,V] = svd(X,'econ');
3+
out = U*shrink(S,tau)*V';
4+
end

‎UTILS/colors2.mat‎

363 Bytes
Binary file not shown.

‎UTILS/coolcolors.mat‎

487 Bytes
Binary file not shown.

‎UTILS/coolcolorsBW.mat‎

701 Bytes
Binary file not shown.

‎UTILS/cosamp.m‎

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,49 @@
1+
function Sest = cosamp(Phi,u,K,tol,maxiterations)
2+
3+
% Cosamp algorithm
4+
% Input
5+
% K : sparsity of Sest
6+
% Phi : measurement matrix
7+
% u: measured vector
8+
% tol : tolerance for approximation between successive solutions.
9+
% Output
10+
% Sest: Solution found by the algorithm
11+
%
12+
% Algorithm as described in "CoSaMP: Iterative signal recovery from
13+
% incomplete and inaccurate samples" by Deanna Needell and Joel Tropp.
14+
%
15+
16+
17+
% This implementation was written by David Mary,
18+
% but modified 20110707 by Bob L. Sturm to make it much clearer,
19+
% and corrected multiple times again and again.
20+
% To begin with, see: http://media.aau.dk/null_space_pursuits/2011/07/ ...
21+
% algorithm-power-hour-compressive-sampling-matching-pursuit-cosamp.html
22+
%
23+
% This script/program is released under the Commons Creative Licence
24+
% with Attribution Non-commercial Share Alike (by-nc-sa)
25+
% http://creativecommons.org/licenses/by-nc-sa/3.0/
26+
% Short Disclaimer: this script is for educational purpose only.
27+
% Longer Disclaimer see http://igorcarron.googlepages.com/disclaimer
28+
29+
% Initialization
30+
Sest = zeros(size(Phi,2),1);
31+
v = u;
32+
t = 1;
33+
numericalprecision = 1e-12;
34+
T = [];
35+
while (t <= maxiterations) && (norm(v)/norm(u) > tol)
36+
y = abs(Phi'*v);
37+
[vals,z] = sort(y,'descend');
38+
Omega = find(y >= vals(2*K) & y > numericalprecision);
39+
T = union(Omega,T);
40+
b = pinv(Phi(:,T))*u;
41+
[vals,z] = sort(abs(b),'descend');
42+
Kgoodindices = (abs(b) >= vals(K) & abs(b) > numericalprecision);
43+
T = T(Kgoodindices);
44+
Sest = zeros(size(Phi,2),1);
45+
b = b(Kgoodindices);
46+
Sest(T) = b;
47+
v = u - Phi(:,T)*b;
48+
t = t+1;
49+
end

‎UTILS/istft/Resynthesis.m‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
clear, clc, close all
2+
3+
% music program (stochastic non-stationary signal)
4+
[x, fs] = wavread('track.wav');
5+
x = x(:, 1);
6+
xmax = max(abs(x));
7+
x = x/xmax;
8+
9+
% signal parameters
10+
xlen = length(x);
11+
t = (0:xlen-1)/fs;
12+
13+
% define analysis and synthesis parameters
14+
wlen = 1024;
15+
h = wlen/4;
16+
nfft = wlen;
17+
18+
% perform time-frequency analysis and resynthesis of the original signal
19+
[stft, f, t_stft] = stft(x, wlen, h, nfft, fs);
20+
[x_istft, t_istft] = istft(stft, h, nfft, fs);
21+
22+
% plot the original signal
23+
figure(1)
24+
plot(t, x, 'b')
25+
grid on
26+
xlim([0 max(t)])
27+
ylim([-1.1 1.1])
28+
set(gca, 'FontName', 'Times New Roman', 'FontSize', 14)
29+
xlabel('Time, s')
30+
ylabel('Normalized amplitude')
31+
title('Original and reconstructed signal')
32+
33+
% plot the resynthesized signal
34+
hold on
35+
plot(t_istft, x_istft, '-.r')
36+
legend('Original signal', 'Reconstructed signal')

‎UTILS/istft/WindowChoice.m‎

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,19 @@
1+
clear, clc, close all
2+
3+
wlen = 64; % window length (recomended to be power of 2), choose it!
4+
h = wlen/4; % hop size (recomended to be power of 2), choose it!
5+
k = 5*wlen; % overlap-add span
6+
win = hamming(wlen, 'periodic'); % window, choose it!
7+
s = zeros(k, 1);
8+
9+
for k = 0:h:k-wlen
10+
indx = k+1:k+wlen; % current window location
11+
s(indx) = s(indx) + win.^2; % window overlap-add
12+
winp(indx) = win; % for plot only
13+
plot(winp, 'ok') % plot just this window
14+
hold on
15+
end
16+
17+
W0 = sum(win.^2); % find W0
18+
s = h*s/W0; % scale the window overlap-add
19+
stem(s, 'r'); % plot window overlap-add

‎UTILS/istft/istft.m‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2+
% Inverse Short-Time Fourier Transform %
3+
% with MATLAB Implementation %
4+
% %
5+
% Author: M.Sc. Eng. Hristo Zhivomirov 12/26/13 %
6+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
7+
8+
function [x, t] = istft(stft, h, nfft, fs)
9+
10+
% function: [x, t] = istft(stft, h, nfft, fs)
11+
% stft - STFT matrix (only unique points, time across columns, freq across rows)
12+
% h - hop size
13+
% nfft - number of FFT points
14+
% fs - sampling frequency, Hz
15+
% x - signal in the time domain
16+
% t - time vector, s
17+
18+
% estimate the length of the signal
19+
coln = size(stft, 2);
20+
xlen = nfft + (coln-1)*h;
21+
x = zeros(1, xlen);
22+
23+
% form a periodic hamming window
24+
win = hamming(nfft, 'periodic');
25+
26+
% perform IFFT and weighted-OLA
27+
if rem(nfft, 2) % odd nfft excludes Nyquist point
28+
for b = 0:h:(h*(coln-1))
29+
% extract FFT points
30+
X = stft(:, 1 + b/h);
31+
X = [X; conj(X(end:-1:2))];
32+
33+
% IFFT
34+
xprim = real(ifft(X));
35+
36+
% weighted-OLA
37+
x((b+1):(b+nfft)) = x((b+1):(b+nfft)) + (xprim.*win)';
38+
end
39+
else % even nfft includes Nyquist point
40+
for b = 0:h:(h*(coln-1))
41+
% extract FFT points
42+
X = stft(:, 1+b/h);
43+
X = [X; conj(X(end-1:-1:2))];
44+
45+
% IFFT
46+
xprim = real(ifft(X));
47+
48+
% weighted-OLA
49+
x((b+1):(b+nfft)) = x((b+1):(b+nfft)) + (xprim.*win)';
50+
end
51+
end
52+
53+
W0 = sum(win.^2); % find W0
54+
x = x.*h/W0; % scale the weighted-OLA
55+
56+
% calculate the time vector
57+
actxlen = length(x); % find actual length of the signal
58+
t = (0:actxlen-1)/fs; % generate time vector
59+
60+
end

‎UTILS/istft/license.txt‎

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,24 @@
1+
Copyright (c) 2015, Hristo Zhivomirov
2+
All rights reserved.
3+
4+
Redistribution and use in source and binary forms, with or without
5+
modification, are permitted provided that the following conditions are
6+
met:
7+
8+
* Redistributions of source code must retain the above copyright
9+
notice, this list of conditions and the following disclaimer.
10+
* Redistributions in binary form must reproduce the above copyright
11+
notice, this list of conditions and the following disclaimer in
12+
the documentation and/or other materials provided with the distribution
13+
14+
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
15+
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16+
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17+
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
18+
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
19+
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
20+
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
21+
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
22+
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
23+
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
24+
POSSIBILITY OF SUCH DAMAGE.

0 commit comments

Comments
 (0)