forked from sivantoledo/javAX25
-
Notifications
You must be signed in to change notification settings - Fork 1
Original MATLAB Experiments
Jose Ricardo de Oliveira Damico edited this page May 30, 2021
·
1 revision
if ~exist('order_multiplier')
order_multiplier = 1
end
% clear java
% javaaddpath e:\files\papers\embedded\ax25
%[x,fs,nbits]=wavread('4Z5LA_3_MAXTRAC.wav');
[x,fs,nbits]=wavread('4Z1PF_3_MAXTRAC.wav');
%s[x,fs,nbits]=wavread('javatest.wav');
%[x,fs,nbits]=wavread('baofeng-new.wav');
%x=xxx;
%fs=48000;
r_fs = [9600, 12000, 11025, 16000, 22050, 24000, 44100, 48000];
r_samples_per_bit = floor(r_fs/1200);
r_nyquist = r_fs/2;
%for r=1:length(r_fs)
% r_b{r}=firls(r_samples_per_bit(r)-1,[0 1050 1150 2250 2350 r_nyquist(r)]/r_nyquist(r) ...
% ,[0 0 1 1 0 0]);
% r_h{r}=fir1(r_samples_per_bit(r)-1,2250/r_nyquist(r));
%end
%filter_lengths = [8 16 32 64 128];
filter_lengths = [1 1];
for l=1:length(filter_lengths)
for r=1:length(r_fs)
f = r_nyquist(r);
%filter_length = filter_lengths(l); % max(r_samples_per_bit(r)-1,1)
filter_length = max(order_multiplier*l*r_samples_per_bit(r)-1,1)
fprintf(1,'// filters of length %d for Fs=%d\n',filter_length,r_fs(r));
%r_b{r,l}=firls(filter_length,[0 1050 1150 2250 2350 f]/f ...
% ,[0 0 1 1 0 0]);
%r_b{r,l}=firls(filter_length,[0 1050 1150 2250 2350 f]/f ...
% ,[0 0 0.33 1 0 0]);
r_b_full{l,r}=firls(filter_length,[0 900 1200 2200 2500 f]/f ...
,[0 0 0.5 1 0 0]);
r_b_none{l,r}=firls(filter_length,[0 900 1200 2200 2500 f]/f ...
,[0 0 1 1 0 0]);
r_b_neg{l,r}=firls(filter_length,[0 900 1200 2200 2500 f]/f ...
,[0 0 1 0.5 0 0]);
r_b_half{l,r}=firls(filter_length,[0 900 1200 2200 2500 f]/f ...
,[0 0 0.75 1 0 0]);
%r_b{r}=firls(r_samples_per_bit(r)-1,[0 800 1000 2400 2600 f]/f ...
% ,[0 0 1 1 0 0]);
%r_b{r}=fir1(filter_length,[ 900 2500]/f,'bandpass');
%r_b{r} = 1;
r_h{l,r}=fir1(filter_length,1200/f);
end
end
r_b = r_b_full;
jf = fopen('e:\files\papers\embedded\javAX25\src\sivantoledo\ax25\Afsk1200Filters.java','w');
if jf==-1
error 'failed to open Java source file'
end
fprintf(jf,'package sivantoledo.ax25;\n');
fprintf(jf,'public class Afsk1200Filters {\n');
fprintf(jf,' static final public int[] sample_rates = {');
for r=1:length(r_fs); fprintf(jf,' %d',r_fs(r)); if r<length(r_fs); fprintf(jf,','); end; end; fprintf(jf,' };\n');
%fprintf(jf,' static final public int[] filter_lengths = {');
%for r=1:length(filter_lengths); fprintf(jf,' %d',filter_lengths(r)); if r<length(filter_lengths); fprintf(jf,','); end; end; fprintf(jf,' };\n');
fprintf(jf,' static final public int[] bit_periods = {');
for r=1:length(r_fs); fprintf(jf,' %d',r_samples_per_bit(r)); if r<length(r_fs); fprintf(jf,','); end; end; fprintf(jf,' };\n');
%=============================
tdf = r_b_full;
fprintf(jf,' static final public float[][][] time_domain_filter_full = {\n');
for l=1:length(filter_lengths)
fprintf(jf,' {\n');
for r=1:length(r_fs);
fprintf(jf,' {');
for k=1:length(tdf{l,r})
fprintf(jf,' %09ef',tdf{l,r}(k));
if k<length(tdf{l,r}); fprintf(jf,','); end;
end
fprintf(jf,' }');
if r<length(r_fs); fprintf(jf,',\n'); else fprintf(jf,'\n'); end;
end;
fprintf(jf,' }');
if l<length(filter_lengths); fprintf(jf,',\n'); else fprintf(jf,'\n'); end
end
fprintf(jf,' };\n');
tdf = r_b_none;
fprintf(jf,' static final public float[][][] time_domain_filter_none = {\n');
for l=1:length(filter_lengths)
fprintf(jf,' {\n');
for r=1:length(r_fs);
fprintf(jf,' {');
for k=1:length(tdf{l,r})
fprintf(jf,' %09ef',tdf{l,r}(k));
if k<length(tdf{l,r}); fprintf(jf,','); end;
end
fprintf(jf,' }');
if r<length(r_fs); fprintf(jf,',\n'); else fprintf(jf,'\n'); end;
end;
fprintf(jf,' }');
if l<length(filter_lengths); fprintf(jf,',\n'); else fprintf(jf,'\n'); end
end
fprintf(jf,' };\n');
%=============================
fprintf(jf,' static final public float[][][] corr_diff_filter = {\n');
for l=1:length(filter_lengths)
fprintf(jf,' {\n');
for r=1:length(r_fs);
fprintf(jf,' {');
for k=1:length(r_h{l,r})
fprintf(jf,' %09ef',r_h{l,r}(k));
if k<length(tdf{l,r}); fprintf(jf,','); end;
end
fprintf(jf,' }');
if r<length(r_fs); fprintf(jf,',\n'); else fprintf(jf,'\n'); end;
end;
fprintf(jf,' }');
if l<length(filter_lengths); fprintf(jf,',\n'); else fprintf(jf,'\n'); end
end
fprintf(jf,' };\n');
fprintf(jf,'}\n');
fclose(jf);
r = 8;
%rate = 9600;
%rate = 11025;
%rate = 8000;
%rate = 48000;
%rate = 24000;
rate = r_fs(r);
%[x,fs,nbits]=wavread('BAOFENG_TR2500_1.wav');
x=x(:,1); % in case it's stereo
decimation = fs/rate
if decimation ~= round(decimation)
%error 'This code can only decimate by integer factors'
%x = decimate(x,decimation);
x = resample(x,rate,fs);
fs = rate;
else
if (decimation ~= 1)
x = x(1:decimation:length(x),1);
fs = fs/decimation;
end
end
if rate==16000
disp('resampling at 8k and interpolating to 16k');
size(x)
x = x(1:2:length(x),1);
% now at 8k
z(1:2:2*length(x),1) = x;
z(2:2:2*length(x)-1,1) = 0.5*(x(1:length(x)-1,1)+x(2:length(x),1) );
%[x,interp_filter] = interp(x,2,4,0.5);
%interp_filter
x = z;
size(x)
end
%plot(z)
%error('xxx')
length(x)
fs
samples_per_bit = floor(fs/1200);
%[samples_per_bit fs/1200]
td_filter_order = samples_per_bit-1;
corrdiff_filter_order = samples_per_bit-1;
%[x,fs,nbits]=wavread('4Z1PF_3_MAXTRAC.wav');
%[x,fs,nbits]=wavread('4Z1PF_3_TR2500.wav');
y=x; % just for testing
nsamples = length(x);
%f1=1050;f2=1150;w1=f1/(fs/2);w2=f2/(fs/2);b=fir1(512,[w1,w2]);
%sy=filter(b,1,x);
%f=900;w=f/(fs/2);highpass=fir1(38,w,'high');
%x = filter(highpass,1,x);
%f1=1150;f2=1250;w1=f1/(fs/2);w2=f2/(fs/2);b1=fir1(samples_per_bit-1,[w1,w2]);
%x=filter(b1,1,x);
%f1=2150;f2=2250;w1=f1/(fs/2);w2=f2/(fs/2);b2=fir1(samples_per_bit-1,[w1,w2]);
%x=filter(b2,1,x);
%f1=1150;f2=2250;w1=f1/(fs/2);w2=f2/(fs/2);b1=fir1(td_filter_order,[w1,w2],'bandpass');
%f1=1250;f2=2150;w1=f1/(fs/2);w2=f2/(fs/2);b2=fir1(td_filter_order,[w1,w2],'stop');
%size(b1)
%size(b2)
%x=filter(b1.*b2,1,x);
%x=filter(b1,1,x);
%nyquist = fs/2;
%b=firls(td_filter_order,[0 1050 1150 2250 2350 nyquist]/nyquist ...
% ,[0 0 1 1 0 0]);
%x=filter(b,1,x);
filter_index = 2;
x = filter(r_b{filter_index,r},1,x);
%f=1300;w=f/(fs/2);h=fir1(corrdiff_filter_order,w);
h = r_h{filter_index,r};
t=0:2*pi/fs:2*pi*nsamples/fs'; t=t(1:nsamples);
t=0:2*pi/fs:2*pi*nsamples/fs'; t=t(1:nsamples);
f0=exp(i*1200*t);
f1=exp(i*2200*t);
%b1
%b2
%h
%f0
%f1
%return;
c0=f0'.*x;
c1=f1'.*x;
z0 = zeros(nsamples,1);
z1 = zeros(nsamples,1);
diff = zeros(nsamples,1);
fdiff = zeros(nsamples,1); % filtered
d = zeros(nsamples,1);
%s = x.^2; %
%s0 = c0.^2;
%s1 = c1.^2;
%periods = nan*zeros(nsamples,1); % symbols energy
%dcd = zeros(nsamples,1); % symbols energy
%gd = zeros(nsamples,1);
%bd = zeros(nsamples,1);
%threshold = 0;
last_transition = 0;
%good_periods_counter = 0;
%bad_periods_counter = 0;
timing_threshold = samples_per_bit / 4;
%bs = []; % bit stream
byte = 0;
bitcount = 0;
%mask=0;
%sample_points = zeros(nsamples,1);
%sample_point = 0;
%stuffcount = 0;
%data = 0;
%databits = 0;
state = 0; % waiting for a packet
%afsk1200 = sivantoledo.ax25.Afsk1200(fs);
%for k=1:length(h)
% fprintf('%09ef,\n',b2(k));
%end
%return
for j=1:nsamples
%afsk1200.addSamples([y(j)],1);
%[f0(j) f1(j) j 11*samples_per_bit]
if j < samples_per_bit || j < length(h)
continue;
end
period = j-samples_per_bit+1:j;
z0(j) = abs(sum(c0(period)));
z1(j) = abs(sum(c1(period)));
diff(j) = z0(j)-z1(j);
fdiff(j) = h*diff(j-length(h)+1:j); % this introduces a delay but from here on we only care about
% zero crossings of fdiff.
% disp(fdiff(j));
d(j) = fdiff(j) > 0;
%dcd(j) = dcd(j-1); % for now; we may flip it
if d(j)~=d(j-1)
p = j - last_transition;
%periods(j) = p;
last_transition = j;
% how many bits?
bits = round(p / samples_per_bit);
%fprintf(' %04f %d %02x %d\n',p/samples_per_bit,bits,byte,bitcount);
if bits==0 || bits>7
state=0; % looking for a flag
continue;
end
%if dcd(j)==1
if bits==7
disp('flag');
data = 0;
databits = 0;
byte = 0;
bitcount = 0;
if true
switch state
case 0
state=1;
case 1
continue;
case 2
%if packet.terminate()
% fprintf(1,'Matlab decoded:\n');
% packet
% fprintf(1,'Matlab decoded end\n');
%else
% disp('bad packet (or no packet)');
%end
disp('done with this packet');
state=1;
otherwise
disp('flag in state other than 0, 1, 2');
end
end
%if state==0; state=1; end; % ready to start a packet
%if state==2; packet.terminate(); state=1; end
else
if true
switch state
case 0
continue; % continue to search for a flag
case 1
state=2; % start a packet
%packet=sivantoledo.ax25.Packet();
case 2
; % fprintf('reading %d bits j=%d\n',bits,j);
end
end
%if state==0 continue;end;
%if state==1; state=2; packet=sivantoledo.ax25.Packet(); disp('starting'); end
for k=1:bits-1
bitcount = bitcount+1;
byte = bitshift(byte,-1);
byte = byte+128; % turn on the least-significant bit
if (bitcount == 8)
fprintf('> %02x %c %c\n',byte,char(byte),char(bitshift(byte,-1)));
%if packet.addByte(byte) == false
% state=0; % looks like a packet, but too long
%end
byte=0;
bitcount=0;
end
end
if bits-1~=5 % the zero we just saw (the transition) was not stuffed
bitcount = bitcount+1;
byte = bitshift(byte,-1);
% byte = byte+1; % turn on the least-significant bit
if (bitcount == 8)
fprintf(1,'> %02x %c %c\n',byte,char(byte),char(bitshift(byte,-1)));
%state
%packet.addByte(byte);
%if packet.addByte(byte) == false
% state=0; % looks like a packet, but too long
%end
byte=0;
bitcount=0;
end
else
fprintf(1,'stuffed bit\n');
end
end
%end % if dcd
end % if d(j) ~= d(j-1)
end % of j loop
clear all;
close all;
clear java
javaaddpath e:\files\papers\embedded\ax25
%======================================================================%
%======================================================================%
order_multiplier = 1;
ax25;
set(0,'defaultfigurecolor','w');
set(0,'defaultfigureInverthardcopy','off');
set(0,'defaultaxescolor','none');
set(0,'defaultfigureinvert','off');
% a bit of very uneven transmission by 4Z5LA
range=72302:73475;
plot(range,y(range),'k-');
a=axis;
axis([72302 73475 a(3) a(4)]);
get(gca)
set(gca,'XTick',[]);
set(gca,'YTick',[]);
set(gca,'Box','off');
set(gca,'XColor',[1 1 1]);
set(gca,'YColor',[1 1 1]);
%line([21398 22240],[0 0],'Color',[0 0 0],'LineStyle','--');
%legend('recorded signal','filtered signal','Location','NorthWest');
print -dtiff -r600 original_uneven.tif
% a segment from 4Z1PF
range=21398:22240;
plot(range,y(range),'k-',range,x(range),'r-');
a=axis;
axis([21398 22240 a(3) a(4)]);
get(gca)
set(gca,'XTick',[]);
set(gca,'YTick',[0]);
line([21398 22240],[0 0],'Color',[0 0 0],'LineStyle','--');
legend('recorded signal','filtered signal','Location','NorthWest');
print -dtiff -r600 original_and_filtered.tif
plot(range,x(range),'r-',range,z0(range),'g-',range,z1(range),'b-');
a=axis;
axis([21398 22240 a(3) a(4)]);
set(gca,'XTick',[]);
set(gca,'YTick',[0]);
line([21398 22240],[0 0],'Color',[0 0 0],'LineStyle','--');
legend('filtered signal','correlation with 1200 Hz','correlation with 2200 Hz','Location','NorthWest');
print -dtiff -r600 filtered_and_correlations.tif
plot(range,x(range),'r-',range,diff(range),'c-',range,fdiff(range),'m-');
a=axis;
axis([21398 22240 a(3) a(4)]);
set(gca,'XTick',[]);
set(gca,'YTick',[0]);
line([21398 22240],[0 0],'Color',[0 0 0],'LineStyle','--');
legend('filtered signal','d = corr(1200)-corr(2200)','d low-pass filtered','Location','SouthWest');
print -dtiff -r600 correlation_diffs.tif
plot(range,fdiff(range),'m-');
a=axis;
mn=min(fdiff(range));
mx=max(fdiff(range));
%axis([21398 22240 a(3) a(4)]);
axis([21398 22240 mn-(mx-mn)*0.3 mx+(mx-mn)*0.15]);
a=axis;
set(gca,'XTick',[]);
set(gca,'YTick',[0]);
line([21398 22240],[0 0],'Color',[0 0 0],'LineStyle','--');
legend('d low-pass filtered','Location','NorthWest');
flips=[];
for t=21398+1:22240
if (fdiff(t-1)*fdiff(t)<=0)
flips=[flips t];
if length(flips)==1
first=t;
end
end
end
%set(gca,'XTick',flips);
%set(gca,'XTick',flips);
ylabelpos = mn-(mx-mn)*0.15;
for i=1:length(flips)
line([flips(i) flips(i)],[a(3) a(4)],'Color',[0 0 0],'LineStyle','--');
if i<length(flips)
periods = (flips(i+1)-flips(i)) / samples_per_bit;
h=text((flips(i+1)+flips(i))/2, ylabelpos ,sprintf('%.3f',periods));
get(h)
set(h,'HorizontalAlignment','center');
set(h,'VerticalAlignment','middle');
set(h,'Rotation',-90);
end
end
print -dtiff -r600 bit_periods.tif
%return
%======================================================================%
%======================================================================%
Fs = 11025;
[H,F]=freqz(r_b_full{2,3},1,2048,Fs);
db=20*log(abs(H))/log(10);
Ht=freqz(r_b_full{2,3},1,[1200 2200],Fs)
Ht=20*log(abs(Ht))/log(10);
desired_x = [0 900 1200 2200 2500 Fs/2];
desired_y = [0 0 0.5 1 0 0];
desired_y = 20*log(abs(desired_y))/log(10);
%plot(F,db,'k-',desired_x,desired_y,'r-');
plot(F,db,'k-');
%get(gca)
set(gca,'YTick',[-120 -80 -40 -6 0]);
set(gca,'XTick',[0 1200 2200 5500]);
axis([min(F) max(F) min(db) 10]);
line([min(F) max(F)],[Ht(1) Ht(1)],'Color',[0 0 0],'LineStyle','--');
line([min(F) max(F)],[Ht(2) Ht(2)],'Color',[0 0 0],'LineStyle','--');
line([1200 1200],[min(db) 10],'Color',[0 0 0],'LineStyle','--');
line([2200 2200],[min(db) 10],'Color',[0 0 0],'LineStyle','--');
title(sprintf('Emphasis and Bandpass Filter of Order %d (Fs=%d)',...
length(r_b_full{2,3}),Fs));
xlabel('Frequency (Hz)');
ylabel('Magnitude Response (dB)');
print -dtiff -r600 emphasis_11025_normal.tif
%======================================================================%
order_multiplier = 10;
ax25;
Fs = 11025;
[H,F]=freqz(r_b_full{2,3},1,2048,Fs);
db=20*log(abs(H))/log(10);
Ht=freqz(r_b_full{2,3},1,[1200 2200],Fs)
Ht=20*log(abs(Ht))/log(10);
plot(F,db,'k-');
%get(gca)
set(gca,'YTick',[-120 -80 -40 -6 0]);
set(gca,'XTick',[0 1200 2200 5500]);
axis([min(F) max(F) min(db) 10]);
line([min(F) max(F)],[Ht(1) Ht(1)],'Color',[0 0 0],'LineStyle','--');
line([min(F) max(F)],[Ht(2) Ht(2)],'Color',[0 0 0],'LineStyle','--');
line([1200 1200],[min(db) 10],'Color',[0 0 0],'LineStyle','--');
line([2200 2200],[min(db) 10],'Color',[0 0 0],'LineStyle','--');
title(sprintf('Emphasis and Bandpass Filter of Order %d (Fs=%d)',...
length(r_b_full{2,3}),Fs));
xlabel('Frequency (Hz)');
ylabel('Magnitude Response (dB)');
print -dtiff -r600 emphasis_11025_highorder.tif
%======================================================================%