function [xd] = recorreder(x,N)
% This function shifts a vector by N positions to the right
% As posted in DSPrelated.com
% http://www.dsprelated.com/showcode/43.php
%
lx = length(x);
if N>lx
fprintf('N has to be smaller than the vector length');
return;
else
xd(N+1:lx) = x(1:lx-N);
end;
clf;
for (truncate = 0:8)
x = imread('aqua.jpg');
subplot(2,2,1);
imshow(x);
subplot(2,2,2);
x_t = (x/2^truncate) * 2^truncate;
%if above line doesn't work your MATLAB is old! use next line:
% x_t = uint8(double(uint8((double(x)/2^truncate))) * 2^truncate);
imshow(x_t);
subplot(2,2,3);
imshow(uint8(abs(double(x_t) - double(x))));
pause(0.5);
end
for zeta = 0.1:0.1:2
sys = tf([9],[1,6*zeta,9]);
p = [1 6*zeta 9];
r = roots(p);
hold on
subplot(1,2,1),plot(real(r), imag(r),'o'), title('Pole Locations'), xlabel('Real'), ylabel('Imaginary')
hold on
subplot(1,2,2),step(sys)
pause(.5)
end
% FFTRESAMPLE Resample a real signal by the ratio p/q
function y = fftResample(x,p,q)
% --- take FFT of signal ---
f = fft(x);
% --- resize in the FFT domain ---
% add/remove the highest frequency components such that len(f) = len2
len1 = length(f);
len2 = round(len1*p/q);
lby2 = 1+len1/2;
if len2 < len1
% remove some high frequency samples
d = len1-len2;
f = f(:);
f(floor(lby2-(d-1)/2:lby2+(d-1)/2)) = [];
elseif len2 > len1
% add some high frequency zero samples
d = len2-len1;
lby2 = floor(lby2);
f = f(:);
f = [f(1:lby2); zeros(d,1); f(lby2+1:end)];
end
% --- take the IFFT ---
% odd number of sample removal/addition may make the FFT of a real signal
% asymmetric and artificially introduce imaginary components - we take the
% real value of the IFFT to remove these, at the cost of not being able to
% reample complex signals
y = real(ifft(f));
function op = vsola(ip, tsm_factor, P)
% Implementation of the variable parameter synchronised overlap add VSOLA algorithm
% This implementation makes use of the standard SOLA algorithm (Rocus and Wilgus, ICASSP 1986) and some efficient paramter settings for the SOLA algorithm (Dorran, Lawlor and Coyle, ICASSP 2003) and (Dorran, Lawlor and Coyle, DAFX 2003)
% Given an input, ip, the longest likely pitch period, P, and and a time scale modification factor, tsm_factor, return an output, op, that is a time scaled version of the input. The synthesis_overlap_size is the length, in samples of the lowest likely fundametal period of the input.
% for speech synthesis_overlap_size is set to (16/1000)*fs samples and for music synthesis_overlap_size is typically set to (20/1000)*fs samples
%
% As with all time-domain algorithms this will only work well for signals which have a strong periodic element. For more complex audio you should use a frequency-domain implementation like the phase vocoder
%
% David Dorran, Audio Research Group, Dublin Institute of Technology
% david.dorran@dit.ie
% http://eleceng.dit.ie/dorran
% http://eleceng.dit.ie/arg
%
% make sure input is mono and transpose if necessary
[r, c] = size(ip);
if r > 1
ip = ip';
end;
[r, c] = size(ip);
if r > 1
disp('Note :only works on mono signals');
op = [];
return
end;
% initialize the values of analysis_frame_length, analysis_window_offset, synthesis_frame_offset and length_of_overlap
desired_tsm_len = round(length(ip)*tsm_factor);
P = round(P); %longest likely pitch period in samples
Lmax = round(P * 1.5);% found this to a reasonable value for the Lmax- Lmax is the duration over which the correlation function is applied
stationary_length = (P * 1.5); % found this to a reasonable value for the stationary length - This is the max duration that could be discarded/repeated
analysis_window_offset = round((stationary_length - P)/abs(tsm_factor - 1)); % this equation was derived.
synthesis_window_offset = round(analysis_window_offset * tsm_factor);
analysis_frame_length = round(Lmax + synthesis_window_offset);
number_of_analysis_frames = floor((length(ip)- analysis_frame_length)/analysis_window_offset);
if number_of_analysis_frames < 2 %not much time-scaling being done just return the input
op = ip;
return;
end;
%the next two lines just ensure that the last frame finishes at the very end of the signal (not essential)
zpad = zeros(1, (number_of_analysis_frames*analysis_window_offset) + analysis_frame_length - length(ip));
ip = [ip zpad];
%initialize the output
op = zeros(1, desired_tsm_len);
%initialize the first output frame
op(1 : analysis_frame_length) = ip(1 : analysis_frame_length);
min_overlap = round(Lmax - P); %ensure that there is some minimum overlap
count = 0;
% Loop for the 2nd analysis frame to the number_of_analysis_frames
for m = 1 : number_of_analysis_frames
%grab the mth input frame
ip_frame = ip(analysis_window_offset * m : (analysis_window_offset * m) + analysis_frame_length - 1);
%grab the maximum overlapping segments from the inout frame and the current output
seg_1 = op(round(synthesis_window_offset*(m-1))+analysis_frame_length - Lmax : round(synthesis_window_offset*(m-1))+analysis_frame_length -1);
seg_2 = ip_frame(1: Lmax);
%compute the correlation of these segments
correlation = xcorr(seg_2, seg_1,'coeff');
%Find the best point to overlap (opt_overlap_length) making sure not to exceed the maximum or go below the minimum overlap.
correlation(length(correlation) - Lmax -1: length(correlation)) = -100;
correlation(1: min_overlap) = -100;
[max_correlation, opt_overlap_length] = max(correlation);
if(max_correlation == 0)
opt_overlap_length = Lmax;
end;
% offset = Lmax - opt_overlap_length;
% if ((offset + analysis_window_offset - synthesis_window_offset) >= 0 & (offset + analysis_window_offset - synthesis_window_offset) <= P)
% count = count +1;
% end;
% append mth analysis frame to the current synthesised output using a linear cross fade
ov_seg = linear_cross_fade(seg_1, ip_frame, opt_overlap_length);
ov_len =(round(synthesis_window_offset*m)+analysis_frame_length) - (round(synthesis_window_offset*(m-1))+analysis_frame_length - Lmax) + 1;
ov_seg = ov_seg(1:ov_len);
op(round(synthesis_window_offset*(m-1))+analysis_frame_length - Lmax: round(synthesis_window_offset*m)+analysis_frame_length) = ov_seg;
end; % end of for loop
% linear cross fade the first segment with the second segment given a certain amount of overlap
% |----------seg1----------|
% |---------------seg2-----------------------|
% |--overlap-|
function op = linear_cross_fade(seg_1, seg_2, overlap_length, cross_fade_duration)
error(nargchk(3,4,nargin));
if nargin < 4
cross_fade_duration = overlap_length;
end
if cross_fade_duration > overlap_length
cross_fade_duration = overlap_length;
end
% overlap the end of seg_1 with the start of seg_2 using a linear cross-fade
if (length(seg_1) < overlap_length),
seg_2 = seg_2(overlap_length - length(seg_1) + 1: length(seg_2));
overlap_length = length(seg_1);
end; % end of if statement
if (length(seg_2) < overlap_length),
seg_1 = seg_1(length(seg_1) - (overlap_length - overlap_length): length(seg_1));
overlap_length = length(seg_2);
end; % end of if statement
overlapping_region = zeros(1, cross_fade_duration); % initialize the overlapping region
seg_1 = seg_1(1: length(seg_1) - (overlap_length - cross_fade_duration));
op = zeros(1, length(seg_1) + length(seg_2) - cross_fade_duration);
if(overlap_length ~= 1)
linear_ramp1 = (cross_fade_duration-1:-1:0) ./(cross_fade_duration-1);
linear_ramp2 = (0:1:cross_fade_duration-1) ./(cross_fade_duration-1);
overlapping_region = (seg_1(length(seg_1)-cross_fade_duration+1:length(seg_1)).*linear_ramp1) + (seg_2(1: cross_fade_duration).*linear_ramp2);
end;
op = [seg_1(1:length(seg_1)-cross_fade_duration) ,overlapping_region , seg_2(cross_fade_duration+1:length(seg_2))];
% END of linear_cross_fade function
[header, s] = hdrload('data_file.txt');
num = 2;
x = zeros(ceil(44100/num),1); %channel 1
y = zeros(ceil(44100/num),1); %channel 2
Fs = 44100; % Sampling frequency
T = 1/Fs; % Sample time
L = Fs/num; % Length of window
g = floor(num*length(s)/Fs); % number of windows
H = 0;
count = 50;
max_index = 1; % pointer for max value
max = 0; % max value
max_index2 = 1; % pointer for max value
max2 = 0; % max value
z = zeros(g,2);
z2 = zeros(g,2);
for i = 1:1:g-3
% creating windows for channel one and two
for j = 1:1:L
k = j + L*i;
x(j,1) = s(k+i,1);
y(j,1) = s(k+i,2);
end
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
YY = fft(y,NFFT)/L;
Y = fft(x,NFFT)/L;
f = Fs/2*linspace(0,1,NFFT/2+1); % creating frequency vector
Y2 = 2*abs(Y(1:NFFT/2+1));
YY2 = 2*abs(YY(1:NFFT/2+1));
%%channel 1
max = Y(max_index,1);
%%channel 2
if(max_index2 < length(YY2))
max2 = YY(max_index2,1);
b = length(YY2) - max_index2 - 2*count;
if b > 0
a = max_index2 + 2*count;
else
a = length(YY2);
end
for j = max_index2:1:(a-1)
if max2 < YY2(j+1,1)
max2 = YY2(j+1,1);
max = Y2(j+1,1);
max_index2 = j+1;
end
end
end
z2(i+1,1) = abs(max2);
z2(i+1,2) = f(1,max_index2);
z(i+1,1) = abs(max);
z(i+1,2) = f(1,max_index2);
if(max_index2 > 4*count)
max_index2 = max_index2 - count;
end
% RECORD MOVIE
% semilogy(z2(:,2),z2(:,1),f,YY2,z(:,2),z(:,1),f,Y2);
% axis([0 20000 10^-3 Inf]);
end
semilogy(z2(:,2),z2(:,1),z(:,2),z(:,1));
axis([0 20000 0.00001 Inf]);
%
% Copyright (C) 2004 by Hernán Ordiales
% <audiocode@uint8.com.ar>
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
function fast_conv( file, impulse, output )
% Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.
% USE: fast_conv( 'source.wav', 'impulse_response.wav', 'output.wav' );
h_stereo = true;
x_stereo = true;
clip_factor = 1.01; % clip factor default value
% impulse response:
[ h, Fs1 ] = wavread( impulse ); % 2 columns matrix, one per channel
h1 = h(:,1); % col 1 -> left channel
h_channels = size(h); h_channels = h_channels(2);
if( h_channels == 2 )
h2 = h(:,2); % col 2 -> right channel
else
h2 = h(:,1); % col 1 -> left channel
h_stereo = false;
end
% file to process:
[ x, Fs2 ] = wavread( file );
x1 = x(:,1); % col 1 -> left channel
x_channels = size(x); x_channels = x_channels(2);
if( x_channels == 2 )
x2 = x(:,2); % col 2 -> right channel
else
x2 = x(:,1); % col 1 -> left channel
x_stereo = false;
end
% if sample rates are the same
if( Fs1 == Fs2 )
% searches for the amount of points required to perform the FFT
L = length(h1) + length(x1) - 1; % linear convolution length
disp('Processing...');
% Use nextpow if you have it to replace below code -> N = 2^nextpow(L);
N = 2;
while( N<L ) % first 2 pow greater than L
N = N*2;
end
% Note: N>=L is needed because the IDFT of the multiplication is the circular convolution and to match it to the
% common one, N>=L is required (where L=N1+N2-1;N1=length(x);N2=length(h))
% FFT(X,N) is the N points FFT, with zero padding if X has less than N points and truncated if has more.
H1 = fft( h1, N ); % Fourier transform of the impulse
X1 = fft( x1, N ); % Fourier transform of the input signal
F_y1 = H1 .* X1; % spectral multiplication
y1 = ifft( F_y1 ); % time domain again
% audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped"
y1 = y1/( max(y1)*clip_factor ); % to avoid clipping
% if one of the files is stereo, process the other channel
if( h_stereo || x_stereo )
H2 = fft( h2, N ); % Fourier transform of the impulse
X2 = fft( x2, N ); % Fourier transform of the input signal
F_y2 = H2 .* X2; % spectral multiplication
y2 = ifft( F_y2 ); % time domain again
y2 = y2/( max(y2)*clip_factor ); % to avoid clipping
y = [ y1 y2 ]; % stereo output (two channels)
else
y = y1; % mono output
end
wavwrite( y, Fs1, output );
disp( 'Convolution success.' );
else
disp('Error: files has different sample rate.');
end
function op = vsola(ip, tsm_factor, P)
% Implementation of the variable parameter synchronised overlap add VSOLA algorithm
% This implementation makes use of the standard SOLA algorithm (Rocus and Wilgus, ICASSP 1986) and some efficient paramter settings for the SOLA algorithm (Dorran, Lawlor and Coyle, ICASSP 2003) and (Dorran, Lawlor and Coyle, DAFX 2003)
% Given an input, ip, the longest likely pitch period, P, and and a time scale modification factor, tsm_factor, return an output, op, that is a time scaled version of the input. The synthesis_overlap_size is the length, in samples of the lowest likely fundametal period of the input.
% for speech synthesis_overlap_size is set to (16/1000)*fs samples and for music synthesis_overlap_size is typically set to (20/1000)*fs samples
%
% As with all time-domain algorithms this will only work well for signals which have a strong periodic element. For more complex audio you should use a frequency-domain implementation like the phase vocoder
%
% David Dorran, Audio Research Group, Dublin Institute of Technology
% david.dorran@dit.ie
% http://eleceng.dit.ie/dorran
% http://eleceng.dit.ie/arg
%
% make sure input is mono and transpose if necessary
[r, c] = size(ip);
if r > 1
ip = ip';
end;
[r, c] = size(ip);
if r > 1
disp('Note :only works on mono signals');
op = [];
return
end;
% initialize the values of analysis_frame_length, analysis_window_offset, synthesis_frame_offset and length_of_overlap
desired_tsm_len = round(length(ip)*tsm_factor);
P = round(P); %longest likely pitch period in samples
Lmax = round(P * 1.5);% found this to a reasonable value for the Lmax- Lmax is the duration over which the correlation function is applied
stationary_length = (P * 1.5); % found this to a reasonable value for the stationary length - This is the max duration that could be discarded/repeated
analysis_window_offset = round((stationary_length - P)/abs(tsm_factor - 1)); % this equation was derived.
synthesis_window_offset = round(analysis_window_offset * tsm_factor);
analysis_frame_length = round(Lmax + synthesis_window_offset);
number_of_analysis_frames = floor((length(ip)- analysis_frame_length)/analysis_window_offset);
if number_of_analysis_frames < 2 %not much time-scaling being done just return the input
op = ip;
return;
end;
%the next two lines just ensure that the last frame finishes at the very end of the signal (not essential)
zpad = zeros(1, (number_of_analysis_frames*analysis_window_offset) + analysis_frame_length - length(ip));
ip = [ip zpad];
%initialize the output
op = zeros(1, desired_tsm_len);
%initialize the first output frame
op(1 : analysis_frame_length) = ip(1 : analysis_frame_length);
min_overlap = round(Lmax - P); %ensure that there is some minimum overlap
count = 0;
% Loop for the 2nd analysis frame to the number_of_analysis_frames
for m = 1 : number_of_analysis_frames
%grab the mth input frame
ip_frame = ip(analysis_window_offset * m : (analysis_window_offset * m) + analysis_frame_length - 1);
%grab the maximum overlapping segments from the inout frame and the current output
seg_1 = op(round(synthesis_window_offset*(m-1))+analysis_frame_length - Lmax : round(synthesis_window_offset*(m-1))+analysis_frame_length -1);
seg_2 = ip_frame(1: Lmax);
%compute the correlation of these segments
correlation = xcorr(seg_2, seg_1,'coeff');
%Find the best point to overlap (opt_overlap_length) making sure not to exceed the maximum or go below the minimum overlap.
correlation(length(correlation) - Lmax -1: length(correlation)) = -100;
correlation(1: min_overlap) = -100;
[max_correlation, opt_overlap_length] = max(correlation);
if(max_correlation == 0)
opt_overlap_length = Lmax;
end;
% offset = Lmax - opt_overlap_length;
% if ((offset + analysis_window_offset - synthesis_window_offset) >= 0 & (offset + analysis_window_offset - synthesis_window_offset) <= P)
% count = count +1;
% end;
% append mth analysis frame to the current synthesised output using a linear cross fade
ov_seg = linear_cross_fade(seg_1, ip_frame, opt_overlap_length);
ov_len =(round(synthesis_window_offset*m)+analysis_frame_length) - (round(synthesis_window_offset*(m-1))+analysis_frame_length - Lmax) + 1;
ov_seg = ov_seg(1:ov_len);
op(round(synthesis_window_offset*(m-1))+analysis_frame_length - Lmax: round(synthesis_window_offset*m)+analysis_frame_length) = ov_seg;
end; % end of for loop
% linear cross fade the first segment with the second segment given a certain amount of overlap
% |----------seg1----------|
% |---------------seg2-----------------------|
% |--overlap-|
function op = linear_cross_fade(seg_1, seg_2, overlap_length, cross_fade_duration)
error(nargchk(3,4,nargin));
if nargin < 4
cross_fade_duration = overlap_length;
end
if cross_fade_duration > overlap_length
cross_fade_duration = overlap_length;
end
% overlap the end of seg_1 with the start of seg_2 using a linear cross-fade
if (length(seg_1) < overlap_length),
seg_2 = seg_2(overlap_length - length(seg_1) + 1: length(seg_2));
overlap_length = length(seg_1);
end; % end of if statement
if (length(seg_2) < overlap_length),
seg_1 = seg_1(length(seg_1) - (overlap_length - overlap_length): length(seg_1));
overlap_length = length(seg_2);
end; % end of if statement
overlapping_region = zeros(1, cross_fade_duration); % initialize the overlapping region
seg_1 = seg_1(1: length(seg_1) - (overlap_length - cross_fade_duration));
op = zeros(1, length(seg_1) + length(seg_2) - cross_fade_duration);
if(overlap_length ~= 1)
linear_ramp1 = (cross_fade_duration-1:-1:0) ./(cross_fade_duration-1);
linear_ramp2 = (0:1:cross_fade_duration-1) ./(cross_fade_duration-1);
overlapping_region = (seg_1(length(seg_1)-cross_fade_duration+1:length(seg_1)).*linear_ramp1) + (seg_2(1: cross_fade_duration).*linear_ramp2);
end;
op = [seg_1(1:length(seg_1)-cross_fade_duration) ,overlapping_region , seg_2(cross_fade_duration+1:length(seg_2))];
% END of linear_cross_fade function
function BaryExample
%Author: J. Selva. 2010.
%
%This function compares the performance of the bary-shifted coefficients with that of the ...
%optimal equirriple coefficients. For a complete discussion see
% J. Selva, 'Design of barycentric interpolator for uniform and nonuniform sampling
% grids', IEEE Trans. on Signal Processing, vol. 58, n. 3, pp. 1618-1627, March 2010.
%This code was used to generate the example in Sec. IV.A of this paper.
T = 1; %Sampling period.
B = 0.7; %Two-sided bandwidth.
P = 5; %Filter semi-length. The total length is 2*P+1.
Lf = 0:B/400:B/2; %Frequency grid for evaluating the maximum ripple.
%Reference time shift, normalized to T.
unRef = 0.5;
%Reference filter coefficients. They introduce a fractional shift 0.5*T.
cRef = cremez(2*P,2*[0,B/2],{@fresp,T,P,T*unRef});
%Instants corresponding to each sample, relative to the central sample.
tg = fliplr((-P:P)*T);
%Grid of normalized (to T) time shifts in which the maximum ripple will be evaluated.
Lun = -P-2:.01:P+2;
%These variables will contain the maximum ripple for each time shift.
LBary = zeros(length(Lun),1);
LRemez = zeros(length(Lun),1);
for k= 1:length(Lun)
LBary(k) = ErrorInFreq(BarycentricShift(cRef,tg,unRef,Lun(k),T),T,P,Lf,Lun(k));
LRemez(k) = ErrorInFreq(cremez(2*P,2*[0,B/2],{@fresp,T,P,T*Lun(k)}),T,P,Lf,Lun(k));
end
%Plot the results.
plot(Lun,mag2db(LRemez),'--',Lun,mag2dB(LBary),'-')
set(gca,'YLim',[-100,0])
xlabel('Normalized time shift (t/T)')
ylabel('Maximum spectral ripple (dB)')
grid on
annotation(gcf,'textarrow',[0.339285714285714 0.251785714285714],...
[0.861904761904766 0.883333333333333],'TextEdgeColor','none',...
'TextBackgroundColor',[1 1 1],...
'String',{'Performance of bary-shifted','coefficients'});
annotation(gcf,'textarrow',[0.521428571428571 0.521428571428571],...
[0.565666666666667 0.495238095238095],'TextEdgeColor','none',...
'TextBackgroundColor',[1 1 1],...
'String',{'For time shifts in [-T,T]','the performance is roughly','the same'});
annotation(gcf,'textarrow',[0.223214285714283 0.214285714285714],...
[0.577571428571438 0.864285714285716],'TextEdgeColor','none',...
'TextBackgroundColor',[1 1 1],...
'String',{'Performance of optimal','equi-ripple coefficients'});
%=========================================================================================
%Applies the barycentric shift.
%cRef is the vector of FIR coefficients that interpolate at instant unRef*T from samples ...
%at instants in tg.
%
%un*T is the desired interpolation instant.
%
%c is the new set of FIR coefficients.
function c = BarycentricShift(cRef,tg,unRef,un,T)
tRef = unRef * T;
t = un * T;
A = sum(cRef);
cRef = cRef / A;
c = cRef .* (tRef-tg) ./ (t-tg);
c = A * c / sum(c);
%=========================================================================================
%Evaluates the maximum spectral ripple of the FIR interpolator. For large P, it is more
%efficient to employ the FFT.
function E = ErrorInFreq(c,T,P,Lf,un)
E = 0;
for k = 1:length(Lf)
s0 = exp(1i*2*pi*Lf(k)*T*(-P:P));
vRef = exp(1i*2*pi*Lf(k)*T*un);
v = c*fliplr(s0).';
E = max([E, abs(v-vRef)]);
end
%Auxiliary function for cremez.
function [DH,DW] = fresp(N,F,GW,W,T,P,u)
if ischar(N)
DH = 'real';
DW = [];
return
end
DH = exp(1i*2*pi*(u-P*T)*GW/2);
DW = ones(size(DH));
%
% Copyright (C) 2004 by Hernán Ordiales
% <audiocode@uint8.com.ar>
%
% This program is free software; you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation; either version 2 of the License, or
% (at your option) any later version.
%
function fast_conv( file, impulse, output )
% Fast convolution of an audio signal with the impulse response of a system modeled as LTI (Linear Time Invariant) one.
% USE: fast_conv( 'source.wav', 'impulse_response.wav', 'output.wav' );
h_stereo = true;
x_stereo = true;
clip_factor = 1.01; % clip factor default value
% impulse response:
[ h, Fs1 ] = wavread( impulse ); % 2 columns matrix, one per channel
h1 = h(:,1); % col 1 -> left channel
h_channels = size(h); h_channels = h_channels(2);
if( h_channels == 2 )
h2 = h(:,2); % col 2 -> right channel
else
h2 = h(:,1); % col 1 -> left channel
h_stereo = false;
end
% file to process:
[ x, Fs2 ] = wavread( file );
x1 = x(:,1); % col 1 -> left channel
x_channels = size(x); x_channels = x_channels(2);
if( x_channels == 2 )
x2 = x(:,2); % col 2 -> right channel
else
x2 = x(:,1); % col 1 -> left channel
x_stereo = false;
end
% if sample rates are the same
if( Fs1 == Fs2 )
% searches for the amount of points required to perform the FFT
L = length(h1) + length(x1) - 1; % linear convolution length
disp('Processing...');
% Use nextpow if you have it to replace below code -> N = 2^nextpow(L);
N = 2;
while( N<L ) % first 2 pow greater than L
N = N*2;
end
% Note: N>=L is needed because the IDFT of the multiplication is the circular convolution and to match it to the
% common one, N>=L is required (where L=N1+N2-1;N1=length(x);N2=length(h))
% FFT(X,N) is the N points FFT, with zero padding if X has less than N points and truncated if has more.
H1 = fft( h1, N ); % Fourier transform of the impulse
X1 = fft( x1, N ); % Fourier transform of the input signal
F_y1 = H1 .* X1; % spectral multiplication
y1 = ifft( F_y1 ); % time domain again
% audio normalization: if "y = y/max(y)" -> "values out of [-1,+1] range are clipped"
y1 = y1/( max(y1)*clip_factor ); % to avoid clipping
% if one of the files is stereo, process the other channel
if( h_stereo || x_stereo )
H2 = fft( h2, N ); % Fourier transform of the impulse
X2 = fft( x2, N ); % Fourier transform of the input signal
F_y2 = H2 .* X2; % spectral multiplication
y2 = ifft( F_y2 ); % time domain again
y2 = y2/( max(y2)*clip_factor ); % to avoid clipping
y = [ y1 y2 ]; % stereo output (two channels)
else
y = y1; % mono output
end
wavwrite( y, Fs1, output );
disp( 'Convolution success.' );
else
disp('Error: files has different sample rate.');
end
clc, clear all, close all
N=128;
x=-1:2/128:1-2/128;
%number of stages
M=q_log2(N);
%to array elements in a proper way for the butterfly
for k=0:N-1;
x1(k+1)=x(bi2de(fliplr(de2bi(k,M)))+1);
end
%to set up twiddle factors
for k=0:N/2-1;
W(k+1)=exp(-1i*2*pi*k/N);
end
%entering to graph
%k holds the stage
for k=1:M;
k2=2^k;
k21=2^(k-1);
Mk=2^(M-k);
%l holds the butterfly
for l=1:N/k2;
%l=1:N/pow(2,k)
%m holds the upper element of the butterfly
for m=1:k21;
op1 = x1((l-1)*(k2)+m);
op2 = x1((l-1)*(k2)+m+k21);
%butterfly equation
x1((l-1)*(k2)+m) = op1+W(Mk*(m-1)+1)*op2;
x1((l-1)*(k2)+m+k21) = op1-W(Mk*(m-1)+1)*op2;
end
end
end
NZ = 1; % number of ZEROS in the filter to be designed
NP = 4; % number of POLES in the filter to be designed
NG = 10; % number of gain measurements
fmin = 100; % lowest measurement frequency (Hz)
fmax = 3000; % highest measurement frequency (Hz)
fs = 10000; % discrete-time sampling rate
Nfft = 512; % FFT size to use
df = (fmax/fmin)^(1/(NG-1)); % uniform log-freq spacing
f = fmin * df .^ (0:NG-1); % measurement frequency axis
% Gain measurements (synthetic example = triangular amp response):
Gdb = 10*[1:NG/2,NG/2:-1:1]/(NG/2); % between 0 and 10 dB gain
% Must decide on a dc value.
% Either use what is known to be true or pick something "maximally
% smooth". Here we do a simple linear extrapolation:
dc_amp = Gdb(1) - f(1)*(Gdb(2)-Gdb(1))/(f(2)-f(1));
% Must also decide on a value at half the sampling rate.
% Use either a realistic estimate or something "maximally smooth".
% Here we do a simple linear extrapolation. While zeroing it
% is appealing, we do not want any zeros on the unit circle here.
Gdb_last_slope = (Gdb(NG) - Gdb(NG-1)) / (f(NG) - f(NG-1));
nyq_amp = Gdb(NG) + Gdb_last_slope * (fs/2 - f(NG));
Gdbe = [dc_amp, Gdb, nyq_amp];
fe = [0,f,fs/2];
NGe = NG+2;
% Resample to a uniform frequency grid, as required by ifft.
% We do this by fitting cubic splines evaluated on the fft grid:
Gdbei = spline(fe,Gdbe); % say `help spline'
fk = fs*[0:Nfft/2]/Nfft; % fft frequency grid (nonneg freqs)
Gdbfk = ppval(Gdbei,fk); % Uniformly resampled amp-resp
figure(1);
semilogx(fk(2:end-1),Gdbfk(2:end-1),'-k'); grid('on');
axis([fmin/2 fmax*2 -3 11]);
hold('on'); semilogx(f,Gdb,'ok');
xlabel('Frequency (Hz)'); ylabel('Magnitude (dB)');
title(['Measured and Extrapolated/Interpolated/Resampled ',...
'Amplitude Response']);
Ns = length(Gdbfk); if Ns~=Nfft/2+1, error("confusion"); end
Sdb = [Gdbfk,Gdbfk(Ns-1:-1:2)]; % install negative-frequencies
S = 10 .^ (Sdb/20); % convert to linear magnitude
s = ifft(S); % desired impulse response
s = real(s); % any imaginary part is quantization noise
tlerr = 100*norm(s(round(0.9*Ns:1.1*Ns)))/norm(s);
disp(sprintf(['Time-limitedness check: Outer 20%% of impulse ' ...
'response is %0.2f %% of total rms'],tlerr));
% = 0.02 percent
if tlerr>1.0 % arbitrarily set 1% as the upper limit allowed
error('Increase Nfft and/or smooth Sdb');
end
figure(2);
plot(s, '-k'); grid('on'); title('Impulse Response');
xlabel('Time (samples)'); ylabel('Amplitude');
c = ifft(Sdb); % compute real cepstrum from log magnitude spectrum
% Check aliasing of cepstrum (in theory there is always some):
caliaserr = 100*norm(c(round(Ns*0.9:Ns*1.1)))/norm(c);
disp(sprintf(['Cepstral time-aliasing check: Outer 20%% of ' ...
'cepstrum holds %0.2f %% of total rms'],caliaserr));
% = 0.09 percent
if caliaserr>1.0 % arbitrary limit
error('Increase Nfft and/or smooth Sdb to shorten cepstrum');
end
% Fold cepstrum to reflect non-min-phase zeros inside unit circle:
% If complex:
% cf = [c(1), c(2:Ns-1)+conj(c(Nfft:-1:Ns+1)), c(Ns), zeros(1,Nfft-Ns)];
cf = [c(1), c(2:Ns-1)+c(Nfft:-1:Ns+1), c(Ns), zeros(1,Nfft-Ns)];
Cf = fft(cf); % = dB_magnitude + j * minimum_phase
Smp = 10 .^ (Cf/20); % minimum-phase spectrum
Smpp = Smp(1:Ns); % nonnegative-frequency portion
wt = 1 ./ (fk+1); % typical weight fn for audio
wk = 2*pi*fk/fs;
[B,A] = invfreqz(Smpp,wk,NZ,NP,wt);
Hh = freqz(B,A,Ns);
figure(3);
plot(fk,db([Smpp(:),Hh(:)])); grid('on');
xlabel('Frequency (Hz)');
ylabel('Magnitude (dB)');
title('Magnitude Frequency Response');
legend('Desired','Filter');
function DemoProposed
% Author: J. Selva.
% First version: Sept. 2007.
%
% This m file contains the code for generating Fig. 2 in Sec. IV.A of paper
% "Convolution-based trigonometric interpolation of band-limited signals", IEEE
% Transactions on Signal Processing, vol. 56, n. 11, pp. 5465-5477, Nov. 2008. The
% equation numbers and references in this m file refer to this publication.
%
% Some comments:
%
% -The function 'KnabAPPeriodicApprox' computes the coefficients G(k/NT) in Sec. IV. This
% function takes several minutes, given that it involves a large number of numerical
% integrals which are performed in function KnabAPSpectrum. However, this must be done
% only once for fixed T, B, P and N, given that the G(k/NT) are independent of the sample
% values.
%
% 'RealTestSignal' implements a band-limited signal which consists of a sum of sinc pulses
% with arbitrary amplitudes and delays.
%
% 'FraserFFTUpsampling' implements the well-known zero-padding FFT
% algorithm. 'ProposedFFTUpsampling' implements the upsampling method in Sec. IV. It is
% worth comparing the code of both functions below.
warning off
T = 1; %Sampling period.
B = 0.5; %Two-sided bandwidth.
P = 18; %Pulse semi-length.
N = 1024; %Number of samples.
%Test signal composed of N+200 sinc pulses with random amplitudes and delays.
Sig = @(t)RealTestSignal(t,B,T,[-100,N+100],N+200);
%Compute coefficients Gn. This takes a while (5 minutes) because it is necessary to
%compute a large number of numerical integrals. However, since this step does not depend
%on the signal samples, it can be performed off-line and once ever.
disp(['Computing spectral coefficients. This takes a few minutes (7 or less), but ', ...
'must be done only',sprintf('\n'),'once ever for fixed T, B, P, and N.'])
% 'KnabAPPeriodicApprox' returns the frequency grid, specified by the first frequency fo
% and the frequency increment Df, the coefficients Gn, one for each grid frequency, and
% the truncation index kg.
[fo,Df,Gn,kg] = KnabAPPeriodicApprox(T,B,P,N);
sRef = Sig((0:N-1)*T); %Vector of signal samples with spacing T.
a = 10; %Oversampling factor.
sUpExact = Sig((0:a*N-1)*T/a); %Exact signal samples with oversampling a.
sUpFraser = FraserFFTUpsampling(sRef,a); %Interpolated samples using a zero-padded FFT.
sUpProposed = ProposedFFTUpsampling(sRef,T,B,12,a,Gn,kg); %Interpolated samples using
%the proposed method.
to = 0; %First instant in time grid.
Dt = T/a; %Grid time increment.
%Plot the interpolation error in dB for the zero-padding FFT method and for the proposed
%method.
plot(to+(0:a*N-1)*Dt,20*log10(abs([sUpExact - sUpFraser;sUpExact - sUpProposed])))
xlabel('Time t/T')
ylabel('Interpolation error (dB)')
annotation(gcf,'textarrow',[0.2911 0.2536],[0.8119 0.7429],...
'String',{'This is the zero-padding','FFT error'});
annotation(gcf,'arrow',[0.3321 0.06964],[0.4833 0.6333]);
annotation(gcf,'textarrow',[0.3446 0.2643],[0.419 0.3548],...
'String',{'This is the error of the','proposed FFT method. Note',...
'that the scale is in dB !!'});
warning on
%=====================================================================
%
% Conventional zero-padding FFT algorithm
function sUp = FraserFFTUpsampling(s,a)
N = length(s);
sF = a * fftshift(fft(s)); %Compute the FFT of the signal samples.
if mod(N,2) == 0 %Pad with zeros to obtain over-sampling factor a.
sF(end) = 0.5*sF(end);
sF = [zeros(1,N*(a-1)/2),sF,zeros(1,N*(a-1)/2)];
else
sF = [zeros(1,(N*(a-1)+1)/2),sF,zeros(1,(N*(a-1)-1)/2)];
end
sUp = ifft(ifftshift(sF)); %Compute the inverse FFT.
%=====================================================================
%
% Proposed FFT algorithm
function sUp = ProposedFFTUpsampling(s,T,B,P,a,Gn,kg)
N = length(s);
sF = fft(s); %Compute the FFT of the signal samples.
sF = sF(1+mod(-kg:kg,N)) .* Gn * (N*a); %Repeat some of the spectral samples periodically,
%and apply the window specified by the vector Gn.
sF = sF(:).';
if mod(N,2) == 0 %Pad with zeros to obtain over-sampling factor a.
sF = [zeros(1,N*a/2-kg),sF,zeros(1,N*a/2-kg-1)];
else
sF = [zeros(1,(N*a-2*kg-1)/2),sF,zeros(1,(N*a-2*kg-1)/2)];
end
sUp = ifft(ifftshift(sF)); %Compute the inverse FFT.
%=====================================================================
%
% Returns the spectral coefficients Gn and the frequency grid specification f0, Df, Gn.
function [f0,Df,Gn,n0] = KnabAPPeriodicApprox(T,B,P,N)
Df = 1/(N*T);
n0 = floor((1/T-B/2)/Df);
f0 = -Df*n0;
Bg = 2/T-B;
Gn = KnabAPSpectrum(f0+Df*(0:2*n0),B,Bg,(P+1/2)*T)/N;
%=====================================================================
%
% Auxiliary function that computes the spectral coefficients Gn.
function G = KnabAPSpectrum(f,B,Bg,TL)
% B: Two-sided signal bandwidth.
% Bg: Two-sided bandwidth of output pulse.
% Be: Two-sided bandwidth of fast-decay pulse.
% TL: Window's width between zeros in time domain.
%
% The two-sided bandwidth of the output pulse is Bg = B+2*Be.
Nf = numel(f);
G = zeros(size(f));
Be = (Bg-B)/2;
Ba = (Bg+B)/2;
Int = @(x)ConvWindow(x,Be,2*TL);
for k = 1:Nf
if f(k)-Ba/2 <= -Be/2 && Be/2 <= f(k)+Ba/2
G(k) = 1;
else
IInt = IntervalIntersection([-Be/2,Be/2],f(k)+[-Ba/2,Ba/2]);
if ~isempty(IInt)
G(k) = quad(Int,IInt(1),IInt(2),1e-16);
end
end
end
% Auxiliary function. It computes samples of the continuous convolution of a square pulse
% with the Kaiser window.
function W = ConvWindow(f,B,TL)
%Kaiser window.
%B: Two-sided window width.
%TL: Width between zeros in time domain.
W = zeros(size(f));
I = abs(f) < B/2;
W(I) = besseli(0,pi*B*(TL/2)*sqrt(1-(2*f(I)/B).^2))/(B*sinc(j*B*TL/2));
function [IInt,RelAInd,RelBInd] = IntervalIntersection(IA,IB)
if IA(2) < IB(1) || IB(2) < IA(1) %No overlap case.
IInt = [];
RelAInd = [];
RelBInd = [];
return
end
SwapIntervals = diff(IA) > diff(IB);
if SwapIntervals %Force IA to be the sorter interval.
Ipr = IA;
IA = IB;
IB = Ipr;
end
if IA(1) < IB(1)
IInt = [IB(1),IA(2)]; %There is intersection to the left.
RelAInd = IB(1)-IA(1);
RelBInd = 0;
if SwapIntervals
Relpr = RelAInd;
RelAInd = RelBInd;
RelBInd = Relpr;
end
return
end
if IA(2) <= IB(2)
IInt = IA; %IA lies inside IB.
RelAInd = 0;
RelBInd = IA(1)-IB(1);
if SwapIntervals
Relpr = RelAInd;
RelAInd = RelBInd;
RelBInd = Relpr;
end
return
end
IInt = [IA(1),IB(2)]; %There is intersection to the right.
RelAInd = 0;
RelBInd = IA(1)-IB(1);
if SwapIntervals
Relpr = RelAInd;
RelAInd = RelBInd;
RelBInd = Relpr;
end
%=====================================================================
%
% Computes samples of a signal formed by sincs with random amplitudes and delays.
function s = RealTestSignal(t,B,T,I,Np)
Sr = rand('seed');
rand('seed',0); %Always select the same random numbers in the sequel.
M = numel(t);
Ap = (rand(Np,1)-0.5)*2 ; %Random amplitudes.
tp = I(1)+rand(Np,1)*diff(I); %Random delays.
rand('seed',Sr)
st = size(t);
t = t(:);
s = zeros(numel(t),1);
for r = 1:Np
s = s + sinc(B*(t-tp(r)))*Ap(r); %Compute sum of sincs
end
s = reshape(s,st);
function plot_amp_spectrum(Fs, y, varargin)
%PLOT_AMP_SPECTRUM Plot amplitude spectrum of a one dimensional signal
% PLOT_AMP_SPECTRUM(FS, Y, COLOR, SEMILOG_SCALE, TITLE_STRING) plots the
% amplitude spectrum of y which has sampling frency Fs. Specification
% of color of the plot, wheather the plot is needed in semilog_scale or
% not, and the title string for the figure are OPTIONAL.
%
% Fs: Sampling frequency.
% y: Plots the one sided amplitude spectrum of the signal with,
% color: color of the plot specified in string,
% semilog_scale: [1]= Plot in log scale.
% Example:
% load ecg;
% plot_amp_spectrum(hz, signal); OR
% plot_amp_spectrum(hz, signal, 'r', 1) OR
% plot_amp_spectrum(hz, signal, 'r', 0, 'Amplitude spectrum of ECG signal')
%--------------------------------------------------------------------------
% Author: Aniket Vartak
% Email: aniket.vartak@gmail.com
% This function uses Matlab's fft() function, and builds on top of it to give
% a handy interface to get a nicely formatted plot of the amplitude
% spectrum.
%--------------------------------------------------------------------------
% Specify the default string incase its not provided
default_title_string = 'Single-Sided Amplitude Spectrum of y(t)';
%Deal with variable number of inputs after fs and x.
if(isempty(varargin)) % No optional information provided
color = 'b';
semilog_scale = 0;
title_string = default_title_string;
end;
if(length(varargin) == 1) % Only color info provided
color = varargin{1};
semilog_scale = 0;
title_string = default_title_string;
end;
if(length(varargin) == 2) % color and integer indicating semilog plot is needed is provided
color = varargin{1};
semilog_scale = varargin{2};
title_string = default_title_string;
end;
if(length(varargin) == 3) % All optional info provided
color = varargin{1};
semilog_scale = varargin{2};
title_string = varargin{3};
end;
%--------------------------------------------------------------------------
T = 1/Fs; % Sample time
L = max(size(y)); % Length of signal
NFFT = 2^nextpow2(L); % Next power of 2 from length of y
Y = fft(y,NFFT)/L; % Take the fourier transform
f = Fs/2*linspace(0,1,NFFT/2+1);% Space out the frequencies
set(gcf, 'Name', strcat('Amplitude Spectrum: [Fs = ',num2str(Fs),' Hz]'));
set(gcf, 'NumberTitle', 'off');
% Plot single-sided amplitude spectrum.
if (semilog_scale)
semilogx(f,20*log10(2*abs(Y(1:NFFT/2+1))), color, 'LineWidth', 2);
else
plot(f,20*log10(2*abs(Y(1:NFFT/2+1))), color,'LineWidth', 2);
end;
title(title_string);
xlabel('Frequency (Hz)');
ylabel('|Y(f)| dB');
grid on;
%==============================================================================
%
% Gray Code Generator for Verilog
%
% This script generates a Verilog header (.vh) file which contains a set of
% of `define directives to establish Gray-coded states for a state machine.
% The output of this script is intended to be used in conjunction with a
% Verilog design which uses Gray coding in some application. For more
% information on Gray coding, please see
% http://mathworld.wolfram.com/GrayCode.html.
%
% Usage:
%
% The bit width of the Gray coded states is set as a script variable.
% 2^bit_width states are generated. The output filename is also
% set as a script variable.
%
% Author: Tom Chatt
%
%------------------------------------------------------------------------------
% Revision History:
%
% Original Version Tom Chatt 11/2/2010
%
%==============================================================================
clear all; clc;
%------------------------------------------------------------------------------
bit_width = 8; % Bit width of the generated code
% Print the generated code entries to a Verilog defines file?
print_to_v_def_file = true;
output_file = ['gray_code_',num2str(bit_width),'.vh']; % Output file name
% Gray code state name prefix. State names will be of the form
% GC<bit_width>B_<state number>.
prefix = ['GC',num2str(bit_width),'B_'];
%--------------Do not edit code below this point-------------------------------
% Generate the code, starting with 1 bit code and looping upward
code = [0; 1];
if bit_width ~= 1
for i = 2 : bit_width
code = [zeros(size(code,1), 1), code; ones(size(code,1), 1), flipud(code)]; %#ok<AGROW>
end
end
% Print out the generated code
fid = fopen(output_file, 'w');
for i = 1 : size(code,1)
macro_name = [prefix, num2str(i)];
macro_val = [num2str(bit_width),'''', dec2bin(binvec2dec(code(i,:)),bit_width)];
fprintf(fid, ['ifndef ',macro_name,'\n']);
fprintf(fid, [' ','`define ', macro_name, ' ', macro_val, '\n']);
fprintf(fid, 'endif\n');
end
fclose('all');
% ----------------------------------------------------------
% Title: Discrete Wavelet Transform
% Author: David Valencia
% UPIITA IPN 2010
%
% For DSPRelated.com
% in: http://www.dsprelated.com/showcode/13.php
%
% Computes the Discrete Wavelet Transform
% equivalent filters
%
% Dependencies:
% formfilter.m - http://www.dsprelated.com/showcode/12.php
% upsample2.m - http://www.dsprelated.com/showcode/10.php
%
% Revisions:
% v1.0a Commented and translated in English
%
% More information at:
% http://www.dsprelated.com/showarticle/115.php
% -----------------------------------------------------------
close all; clear all; clc;
disp('== GENERALIZED DWT FILTERS ==')
%% STEP 1 - Base Filters
% Select type of base filters
typeofbasis = 'o';
typbior = 'bior2.2';
typor = 'db6';
% Obtain base filters
if(typeofbasis == 'b')
[Rf,Df] = biorwavf(typbior);
[h0,h1,g0,g1] = biorfilt(Df,Rf);
elseif (typeofbasis == 'o')
[h0,h1,g0,g1] = wfilters(typor);
end;
%% STEP 2 - Parameter configuration
% One can declare or recover an input vector from another
% program, but here an example vector is built for testing.
L = length(h0); %Base filter lenght
n_stages = 2; %Of the low-pass stage
n_branches = 2^n_stages; %Number of branches
% L = Basic filters length (h0 ó h1)
% N = Input vector length
% n_stages = stage quantity, it generates (2^n_stages) branches
figure;
[H1,Wg]=freqz(h1,1,512);
plot(Wg,abs(H1),'r','linewidth',2); hold on;
for i = 0:(n_branches/2)-1
if n_stages < 3
hx = formfilter(n_stages, (n_branches/2) - i, h0, h1);
else
hx = formfilter(n_stages, (n_branches/2) - i - 1, h0, h1);
end
Lhx = length(hx);
% Graphing code
[H0,Wg]=freqz(hx,1,512);
plot(Wg,abs(H0),'b','linewidth',2); hold on;
pause; % Used to see how each filter appears
end
