Adaptive Filter | ADPCM | ADSP | ADSP-2181 | Aliasing | AMR | Anti-Aliasing | ARMA | Autocorrelation | AutoCovariance | Beamforming | Bessel | Blackfin | Butterworth | C6713 | CCS | Chebyshev | CIC Filter | Circular Convolution | Code Composer Studio | Comb Filter | Compression | Convolution | Cross Correlation | DCT | Decimation | Deconvolution | Demodulation | DM642 | DSP Boards | DSP/BIOS | DTMF | Echo Cancellation | Equalization | Equalizer | ETSI | EZLITE (Ez-kit Lite) | FFT | FFTW | FIR Filter | Fixed Point | FSK | G.711 | G.723 | G.729 | Gaussian Noise | Goertzel | GPIO | Hilbert Transform | IFFT | IIR Filter | Interpolation | Invariance | JTAG | Kalman | Laplace Transform | Levinson | LPC | McBSP | MIPS | Modulation | MPEG | Multirate | Notch Filter | Nyquist | OFDM | Oversampling | Pink Noise | Pitch | PLL | Polyphase | QAM | QDMA | Quantization | Quantizer | Radar | Random Noise | Reed Solomon | Remez | Resampling | RTDX | Sampling | Sharc | TI C6711 | Undersampling | Viterbi | Wavelets | White Noise | Wiener Filter | Windowing | XDS510PP | Z Transform

Hi,

I'm interested in warped filtering and have had a look at some tutorials.

I tried to implement a little program that warps a sine signal along the
frequency axis. Unfortunately, it doesn't seem to work as intended (see code
below)

Does anyone have a (efficient) way of implementing a warped delay chain?

Cheers
Thomas


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
clear all;clc;close all;

% Define time vector - fs = 44.1e3 kHz
t = 0:1/44.1e3:(1000*1/44.1e3-1/44.1e3);
% Time signal to be frequency-warped - 500 Hz sine
x1 = sin(2*pi*500*t);

% Calculate power spectrum
[Pxx,F]=psd(x1,length(t),441e3,hanning(length(t)));

% plot inputs in time and frequency
figure(1),
subplot(2,1,1),plot(t,x1);set(gca,'xlim',[min(t) max(t)]);grid;title('500
Hz');
subplot(2,1,2),plot(F,Pxx);grid;

% Plot group delay for a single delay element
figure(2),
a    = [0.7],
num  = [-a 1];
den  = fliplr(num);
[G,W]= grpdelay(num,den,512,44.1e3);
plot(W,G);title('group-delay');grid;

% filter time signal with chain of delay elements
y1 = x1;
no_iter = 1000;

for indx=1:no_iter,
    y1 = filter(num,den,y1);
    out1(:,indx)=y1';
end
[Pyy,F]=psd(out1(:,1),256,44.1e3,hanning(256))

figure(3),
t = 1:length(y1);
subplot(3,1,1),plot(t,x1,'-k');
subplot(3,1,2),plot(t,out1(:,1),'-k');
subplot(3,1,3),plot(F,Pyy,'-k');

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

I found a few bugs in my matlab program.....

At glance it seems to frequency-warp the sine signal but it looses
amplitude - any suggestions?


clear all;clc;close all;

% Define time vector
t = 0:1/44.1e3:(1000*1/44.1e3-1/44.1e3);
% time signal to be frequency-warped
x1 = sin(2*pi*500*t);

% Calculate power spectrum
[Pxx,F]=psd(x1,length(t),441e3,hanning(length(t)));

% plot inputs in time and frequency
figure(1),
subplot(2,1,1),plot(t,x1);set(gca,'xlim',[min(t) max(t)]);grid;title('500
Hz');
subplot(2,1,2),plot(F,Pxx);grid;

% Plot group delay for a single delay element
figure(2),
a    = [0.1],
num  = [-a 1];
den  = fliplr(num);
[G,W]= grpdelay(num,den,512,44.1e3);
plot(W,G);title('group-delay');grid;

% filter time signal with chain of delay elements
y1 = x1;
no_iter = 1000;

for indx=1:no_iter,
    y1 = filter(num,den,y1);
    out1(:,indx)=y1';
end
[Pyy,F]=psd(diag(out1),256,44.1e3,hanning(256))

figure(3),
t = 1:length(y1);
subplot(3,1,1),plot(t,x1,'-k');
subplot(3,1,2),plot(t,diag(out1),'-k');
subplot(3,1,3),plot(F,Pyy,'-k');

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

"heureka" <s...@hotmail.com> wrote in message news:<cd1851$1e2d$1...@news.cybercity.dk>...
> I found a few bugs in my matlab program.....
> 
> At glance it seems to frequency-warp the sine signal but it looses
> amplitude - any suggestions?
> 
> clear all;clc;close all;
> 
> % Define time vector
> 
> % filter time signal with chain of delay elements
> y1 = x1;
> no_iter = 1000;
> 
> for indx=1:no_iter,
>     y1 = filter(num,den,y1);
>     out1(:,indx)=y1';
> end
> [Pyy,F]=psd(diag(out1),256,44.1e3,hanning(256))
> 

IMHO, your filter is an all pass filter but still introduces delay
with "filter" command. To see this, check the out1 matrix. So in the
above loop, y1 is essentially the filtered delay version of previous
y1. Since your x1 has 1000 samples and out1 also 1000x1000, so the
delay supresses the magnitude. Can you explain why only takes the
diagonal elements of out1? Because I think if you only take the
diagonal elements of out1, the magnitude should be reducing.

In order to fix it, you may want to design your own "filter" command
which is a convolution operation, I mean, fix the delay after each
iteration in the loop.

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

> IMHO, your filter is an all pass filter but still introduces delay
> with "filter" command. To see this, check the out1 matrix. So in the
> above loop, y1 is essentially the filtered delay version of previous
> y1. Since your x1 has 1000 samples and out1 also 1000x1000, so the
> delay supresses the magnitude. Can you explain why only takes the
> diagonal elements of out1? Because I think if you only take the
> diagonal elements of out1, the magnitude should be reducing.
>
> In order to fix it, you may want to design your own "filter" command
> which is a convolution operation, I mean, fix the delay after each
> iteration in the loop.

Hi Steve,

At first I thought that the diagonal elements were the output of each "delay
tap" in the chain, but that assumption was not correct. The results is at
row number 1000 instead.

Please check out the new script that seem to work. The implementation still
takes quite a lot time, so for the fun of it it could be interesting to see
a faster implementation.

Thanks for your reply
Thomas S. :)

clear all;clc;close all;

% Define time vector
t = 0:1/44.1e3:(1000*1/44.1e3-1/44.1e3);
% time signal to be frequency-warped
x1 = sin(2*pi*500*t);

% Calculate power spectrum
[Pxx,F]=psd(x1,length(t),441e3,hanning(length(t)));

% plot inputs in time and frequency
figure(1),
subplot(2,1,1),plot(t,x1);set(gca,'xlim',[min(t)
max(t)]);grid;title('500Hz');
subplot(2,1,2),plot(F,Pxx);grid;

% Plot group delay for a single delay element
figure(2),
a    = [0.723],
num  = [-a 1];
den  = fliplr(num);
[G,W]= grpdelay(num,den,512,44.1e3);
plot(W,G);title('group-delay');grid;

% filter time signal with chain of delay elements
y1 = x1;
no_iter = 1000;

for indx=1:no_iter,
    y1 = filter(num,den,y1);
    out1(:,indx)=y1';
end
[Pyy,F]=psd(diag(out1),256,44.1e3,hanning(256))

figure(3),
t = 1:length(y1);
subplot(3,1,1),plot(t,x1,'-k');
subplot(3,1,2),plot(t,out1(1000,:),'-k');
subplot(3,1,3),plot(F,Pyy,'-k');

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

No it's still wrong

The envelope is wrong......

for example:

clear all;clc;close all;

% Define time vector
t = 0:1/44.1e3:(1000*1/44.1e3-1/44.1e3);
% time signal to be frequency-warped
x1 = sin(2*pi*500*t).*linspace(0,1,length(t));

% Calculate power spectrum
[Pxx,F]=psd(x1,length(t),441e3,hanning(length(t)));

% plot inputs in time and frequency
figure(1),
subplot(2,1,1),plot(t,x1);set(gca,'xlim',[min(t)
max(t)]);grid;title('500Hz');
subplot(2,1,2),plot(F,Pxx);grid;

% Plot group delay for a single delay element
figure(2),
a    = [0.0723],
num  = [-a 1];
den  = fliplr(num);
[G,W]= grpdelay(num,den,512,44.1e3);
plot(W,G);title('group-delay');grid;

% filter time signal with chain of delay elements
y1 = x1;
no_iter = 500;

for indx=1:no_iter,
    y1 = filter(num,den,y1);
    out1(:,indx)=y1';
end
[Pyy,F]=psd(diag(out1),256,44.1e3,hanning(256))

figure(3),
t = 1:length(y1);
subplot(3,1,1),plot(x1,'-k');
subplot(3,1,2),plot(out1(no_iter,:),'-k');
subplot(3,1,3),plot(F,Pyy,'-k');

figure(4),
subplot(1,2,1),imagesc(out1);grid;
subplot(1,2,2),plot(out1(no_iter,:));grid;

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.

"JohnDoe" <s...@hotmail.com> wrote in message news:<cd49vh$2bvt$1...@news.cybercity.dk>...
> No it's still wrong
> 
> The envelope is wrong......
> 
> for example:
> 
> clear all;clc;close all;
> 
> % Define time vector
> t = 0:1/44.1e3:(1000*1/44.1e3-1/44.1e3);
> % time signal to be frequency-warped
> x1 = sin(2*pi*500*t).*linspace(0,1,length(t));
> 
> % Calculate power spectrum
> [Pxx,F]=psd(x1,length(t),441e3,hanning(length(t)));
> 
> % plot inputs in time and frequency
> figure(1),
> subplot(2,1,1),plot(t,x1);set(gca,'xlim',[min(t)
> max(t)]);grid;title('500Hz');
> subplot(2,1,2),plot(F,Pxx);grid;
> 
> % Plot group delay for a single delay element
> figure(2),
> a    = [0.0723],
> num  = [-a 1];
> den  = fliplr(num);
> [G,W]= grpdelay(num,den,512,44.1e3);
> plot(W,G);title('group-delay');grid;
> 
> % filter time signal with chain of delay elements
> y1 = x1;
> no_iter = 500;
> 
> for indx=1:no_iter,
>     y1 = filter(num,den,y1);
>     out1(:,indx)=y1';
> end
> [Pyy,F]=psd(diag(out1),256,44.1e3,hanning(256))
> 
> figure(3),
> t = 1:length(y1);
> subplot(3,1,1),plot(x1,'-k');
> subplot(3,1,2),plot(out1(no_iter,:),'-k');
> subplot(3,1,3),plot(F,Pyy,'-k');
> 
> figure(4),
> subplot(1,2,1),imagesc(out1);grid;
> subplot(1,2,2),plot(out1(no_iter,:));grid;

If the last row is used, then that makes sense. In this code, what if
you change no_iter to 1000, ie., the same length as the input data,
what do you get?

______________________________
New DSP Code Snippets Section now Live.   Learn more about the reward program for contributors here.