0
$\begingroup$

I am working on signals to check how approximation can change the results. I am testing on example that can change frequency from 22kHz to 20 KHz. I prepare this script in MatLab

clc;
close all;
clear all;
Fs1 = 22000;
Fs2 = 20000;
TimeStep = 1/Fs1;
Length = 100000; 
time = (0:Length-1)*TimeStep;
InputSignal = sin(2*pi*8000*time);
nGcd = gcd(Fs1, Fs2);
wstop = Fs2/(Fs2/nGcd*Fs1);
wpass = 0.9 * wstop;
A = [1 1 0 0];
F = [0 wpass wstop 1];
FilterOrder = 255;
firCoeff = fir2(FilterOrder,F,A);
FilterOrder = length(firCoeff);
BigStep1 =  1000000000/Fs1; %this is 45455,55
BigStep2 = 45455;
Buffer24 = Timestamps(Length, BigStep1); %buffers of timestamps
Buffer24wrong = Timestamps(Length, BigStep2);
SmallerStep = 50000;
Length2 = Length*BigStep2/SmallerStep;
Buffer20 = Timestamps(round(Length2), SmallerStep); %output timestamps
upY = upsample(InputSignal,Fs2/nGcd);
ApplySignal = applyFilter2(FilterOrder,firCoeff,upY,Buffer20,Buffer24);
FilteredSignal= downsample(ApplySignal,Fs1/nGcd);
Length = length(FilteredSignal);
NFFT = 2^nextpow2(Length);
Y = fft(FilteredSignal,NFFT)/Length;
f = Fs2/2*linspace(0,1,NFFT/2+1);
FilteredSignal2 = applyFilter2(FilterOrder,firCoeff,InputSignal,Buffer20,Buffer24wrong);
LengthW = length(FilteredSignal2);
NFFT2 = 2^nextpow2(LengthW);
Y2 = fft(FilteredSignal2,NFFT2)/Length;
f2 = Fs2/2*linspace(0,1,NFFT2/2+1);
figure
subplot(2,2,1)
plot(f,2*abs(Y(1:NFFT/2+1)))
hold on
title('Wrong Input with up&dw sampling')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
subplot(2,2,2)
plot(f2,2*abs(Y2(1:NFFT2/2+1)))
hold on
title('Good inpuit olny without anysampling')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
subplot(2,2,[3,4]);
Index1 = GetIndex(Buffer24,Buffer20(5530),1);
Index2 = GetIndex(Buffer24,Buffer20(5540),1);
plot(Buffer24wrong(Index1:Index2),InputSignal(Index1:Index2),'-o',Buffer24(Index1:Index2),InputSignal(Index1:Index2),'-o', Buffer20(5530:5540),FilteredSignal(5530:5540),'d');
hold on
set(gca,'XTick',Buffer20(5530):SmallerStep:Buffer20(5540))
grid on 
legend('Wrong signal','Up/Dwsampled Signal', 'Output Timestamps')
title('Input Signals')
xlabel('Time')
ylabel('|Y(f)|')

An also function GetIndex:

function Index = GetIndex(EquidistancyBuffer, ReqTime, StartIndex)
LastIndex = length(EquidistancyBuffer);

while((LastIndex-StartIndex)>1)
    Mid = round((StartIndex+LastIndex)/2);
    ActualTime = EquidistancyBuffer(Mid);
    if (ReqTime==ActualTime)
       LastIndex = StartIndex; 
    elseif (ReqTime<ActualTime)
        LastIndex = Mid;
    else
        StartIndex = Mid;
    end

end
Index = Mid;

end

Function applyFilter2:

function yOut = applyFilter2(FilterOrder, Coeff, Signal, ReqTimestmps,EquidistancyBuffer)

StopIndex=length(ReqTimestmps);
for i = FilterOrder+1: StopIndex
    if(i == FilterOrder+1)
        startIndex = 1;
    else
        startIndex = i-1;
    end
    Index = GetIndex(EquidistancyBuffer,ReqTimestmps(i),startIndex);
    yOut(i-FilterOrder) = OutputCalculation(Index-255,FilterOrder,Coeff,Signal);

end

and Timestamps:

function Timestamps_Buffer = Timestamps(Size, Step)

for i = 1: Size
    Timestamps_Buffer(i)=(i-1)*Step;
end
end

I expected on the left subplot only peak on almost 8kHz but I have one closer to 9KHz. When I change Bigstep2 to equal to BigStep1 the results are almost the same. Please explain me what I did wrong. I attached the result of calcultaion also.enter image description here

$\endgroup$
  • $\begingroup$ This question is very close to a programming question. You are effectively asking us to debug your code, which is generally better asked on SO itself. Why don't you just do m=lcm(fsin,fsout); up=m/fsin; down=m/fsout; y=resample(x,up,down); to resample? It seems like you are trying to reinvent the wheel. $\endgroup$ – Peter K. Sep 3 '15 at 15:38
  • $\begingroup$ Oh, come on! By design it has to be upsamled, filtered and downsampled. $\endgroup$ – eryk Sep 4 '15 at 6:25
  • $\begingroup$ All of which is already done in resample. You should read about it. Why bother redoing it? I understand you may want to for pedagogical reasons, but as I said... your question, as it stands, is a programming one not a signal processing one. $\endgroup$ – Peter K. Sep 4 '15 at 12:19
  • $\begingroup$ You are not helping, my friend. It is really about dsp not programming man, the script is OK but the peak is expected on 8kHz not almost 9kHz and this is processing problem. I have the same problem for example with anti aliasing from 25000 Hz to 20000Hz. I have a peak in 5KHz after resampling but expected should be in 8kHz. $\endgroup$ – eryk Sep 4 '15 at 13:15
  • $\begingroup$ Apparently the script is not OK, or you wouldn't be asking the question? I suggest you look at each individual step in your procedure. Those errant frequencies are coming in somewhere in your code. For example, what do InputSignal, upY, and ApplySignal look like? $\endgroup$ – Peter K. Sep 4 '15 at 13:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.