updated pfb.wrap

.gitignore

+bin_dat/
+bin_dat_ack/
 install/
 build/
 subprojects/grrr
 subprojects/pfb
 *.bin
-sadtest.cu
\ No newline at end of file
+sadtest.cu
+*.png
+*.txt
+*.zip
+*.jpg
\ No newline at end of file

compare_input.py

+import struct
+import matplotlib.pyplot as plt 
+from scipy import signal, fft
+import numpy as np
+
+nblocks=2
+read_out_len = 38707200//2 #in_len*nblocks*2
+
+with open("in_6_v2.bin",mode='rb') as file:
+    output_data_f1 = np.fromfile(file, dtype=np.int8, count=read_out_len)
+    #output_data_f1 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+    input_data_c1 = np.array([ complex(output_data_f1[i], output_data_f1[i+1]) for i in range(0,read_out_len,2) ] ,dtype=complex)
+
+
+with open("in_7_v2.bin",mode='rb') as file:
+    # output_data_f2 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+    output_data_f2 = np.fromfile(file, dtype=np.int8, count=read_out_len)
+    input_data_c2 = np.array([ complex(output_data_f2[i], output_data_f2[i+1]) for i in range(0,read_out_len,2) ] ,dtype=complex)
+    
+
+diff = abs(output_data_f2-output_data_f1)
+print(max(diff))
+if max(diff) > 0:
+    idx = np.where(diff==max(diff))
+    idx=idx[0]
+    print(idx)
+    print(output_data_f1[idx])
+    print(output_data_f2[idx])
+    Nc = 25
+    Ni = 144
+    stride = (Nc*Ni)*2
+    print(stride)
+    # plt.plot(idx[0:stride*10],output_data_f1[idx[0:stride*10]],'*')
+    # plt.plot(idx[0:stride*10],output_data_f2[idx[0:stride*10]],'*')
+    # plt.plot(output_data_f1[idx[0]::stride])
+    # plt.plot(output_data_f2[idx[0]::stride])
+    N=100
+    plt.plot(range(idx[0],idx[0]+stride*N,stride),output_data_f1[idx[0]:idx[0]+stride*N:stride])
+    plt.plot(range(idx[0],idx[0]+stride*N,stride),output_data_f2[idx[0]:idx[0]+stride*N:stride])
+    plt.savefig("in.png")
+
+    diff = abs(abs(input_data_c1) -abs(input_data_c2))
+    idx = np.where(diff==max(diff))
+    idx=idx[0]
+    print(idx)
+    print(input_data_c1[idx[0]],input_data_c2[idx[0]])
+
+
+
+# print(input_data_c1[0],input_data_c2[0])
+# print(input_data_c1[0:stride*20:stride])
+# print(input_data_c2[0:stride*20:stride])
+    # stride = nbin*nbeam
+
+    # if plot_on:
+        
+    #     for t in range(0,len(times)):
+
+    #         for c in range(0,1):
+    #             offset = beam_to_plot + nbeam*c + times[t]*nbeam*nbin*nfft + block_to_plot*bf_len
+    #             indexes = np.arange(offset,offset+nfft*stride,stride,dtype=int)
+
+    #             time_in_data_v2 = input_data_c2[indexes]
+    #             print(c, max(abs(time_in_data_v2)))
+    #             print(abs(time_in_data_v2[0:20]))
+
+    #             axs[0][t].plot((abs(time_in_data_v2)))
+    #             axs[0][t].set_xlabel('Frequency bin')
+    #             axs[0][t].set_ylabel('Magnitude')
+    #             axs[0][t].set_title('PFB output t(' + str(times[t]) + ')')
+    #             axs[0][t].grid()
+                
+        

meson_options.txt

@@ -5,7 +5,7 @@ option('RTBF_NBEAM',         type: 'integer', min : 8,   max : 1024,  value : 80
 option('RTBF_STI_LENGTH',    type: 'integer', min : 8,   max : 1024,  value : 21,   description: 'short term integration length')
 option('RTBF_STI_GPU_BLOCK', type: 'integer', min : 8,   max : 1024,  value : 32,   description: 'sti per cuda block')
 option('BIN_WITH_DATA', type: 'integer', min : 0,   max : 25,  value : 1,   description: 'put tone in this channel/bin')
-option('TONE_FREQ', type: 'integer', min : 0,   max : 1024,  value : 1,   description: 'choose the frequency of the tone')
+option('TONE_FREQ', type: 'integer', min : 0,   max : 256,  value : 128,   description: 'choose the frequency of the tone')
 
 
 option('CUDA_ARCH', type:'string', value:'sm_86', description: 'cuda architecture')

pfb_plot.py

@@ -18,13 +18,58 @@ ntime = 5376
 nbeam = 80
 nbin = 25
 fs = 500e6 
-beam_to_plot = 8
+beam_to_plot = 0
 block_to_plot = 0
+bf_len = ntime * nbin * nbeam
+nfft = 256
 
+f, axs = plt.subplots(3, 2,  constrained_layout = True)
+try:
+    print("checking PFB output")
+    read_out_len = bf_len*2*nblocks
+    
+    with open("pfb_6_v2.bin",mode='rb') as file:
+        output_data_f2 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+        input_data_c2 = np.array([ complex(output_data_f2[i], output_data_f2[i+1]) for i in range(0,read_out_len,2) ] ,dtype=complex)
+        stride = nbin*nbeam
+        print(output_data_f2[0],input_data_c2[0])
+        if plot_on:
+            
+            for c in range(0,1):
+                offset = beam_to_plot + nbeam*c + block_to_plot*bf_len
+                indexes = np.arange(offset,offset+nfft*stride,stride,dtype=int)
+
+                time_in_data_v2 = input_data_c2[indexes]
+                print(c, max(abs(time_in_data_v2)))
+
+                axs[0][0].plot(abs(time_in_data_v2))
+                axs[0][0].set_xlabel('Frequency bin')
+                axs[0][0].set_ylabel('Magnitude')
+                axs[0][0].set_title('abs')
+                axs[0][0].grid()
+                # axs[0][0].set_ylim([0.04,0.06])
+
+                axs[1][0].plot(np.real(time_in_data_v2))
+                axs[1][0].set_xlabel('Frequency bin')
+                axs[1][0].set_ylabel('Magnitude')
+                axs[1][0].set_title('real')
+                axs[1][0].grid()
+
+                axs[2][0].plot(np.imag(time_in_data_v2))
+                axs[2][0].set_xlabel('Frequency bin')
+                axs[2][0].set_ylabel('Magnitude')
+                axs[2][0].set_title('imag')
+                axs[2][0].grid()
+
+            
+            f.suptitle('PFB & STI output data')
+
+except Exception as e:
+    print(e)
 
 try:
     print("checking fine STI output")
-    nfft = 256
+    
     out_len_streams = nfft * nbin * nbeam//2 * 4 #4 floats / stivec
     read_out_len = nblocks * out_len_streams
 
@@ -36,79 +81,87 @@ try:
         output_data_c2 = np.array([ complex(output_data_f1[i], output_data_f1[i+1]) for i in range(2,out_len_streams,4) ] )
 
     if plot_on: 
-        f, axs = plt.subplots(3, 1,  constrained_layout = True)
        
         f.set_size_inches(6, 6)
-        for c in range(0,25):
+        for c in range(0,1):
             offset = beam_to_plot + c*nbeam//2
             indexes = np.arange(offset,offset+nfft*stride,stride,dtype=int)
             
             time_in_data = output_data_0[indexes]
-            axs[0].stem(np.arange(0,nfft)+c*nfft, time_in_data, 'C'+str(c % 2)+'o')
-            axs[0].set_title("XX*")
-            axs[0].set_ylabel("magnitude")
-            axs[0].set_xlabel("freq bin")
+            # print(c,time_in_data)
+            axs[0][1].plot(np.arange(0,nfft)+c*nfft, time_in_data)#, 'C'+str(c % 2)+'o')
+            axs[0][1].set_title("XX*")
+            axs[0][1].set_ylabel("magnitude")
+            axs[0][1].set_xlabel("freq bin")
+            # axs[0][1].set_ylim([-0.01,0.01])
+            axs[0][1].grid()
 
             time_in_data = output_data_1[indexes]
-            axs[1].stem(np.arange(0,nfft)+c*nfft, time_in_data, 'C'+str(c % 2)+'o')
-            axs[1].set_title("YY*")
-            axs[1].set_ylabel("magnitude")
-            axs[1].set_xlabel("freq bin")
+            axs[1][1].plot(np.arange(0,nfft)+c*nfft, time_in_data)#, 'C'+str(c % 2)+'o')
+            axs[1][1].set_title("YY*")
+            axs[1][1].set_ylabel("magnitude")
+            axs[1][1].set_xlabel("freq bin")
+            # axs[1][1].set_ylim([-0.01,0.01])
+            axs[1][1].grid()
 
             time_in_data = output_data_c2[indexes]
-            axs[2].stem(np.arange(0,nfft)+c*nfft, np.abs(time_in_data), 'C'+str(c % 2)+'o')
-            axs[2].set_title("XY*+jYX*")
-            axs[2].set_ylabel("magnitude")
-            axs[2].set_xlabel("freq bin")
+            axs[2][1].plot(np.arange(0,nfft)+c*nfft, np.abs(time_in_data))#, 'C'+str(c % 2)+'o')
+            axs[2][1].set_title("XY*+jYX*")
+            axs[2][1].set_ylabel("magnitude")
+            axs[2][1].set_xlabel("freq bin")
+            # axs[2][1].set_ylim([-0.01,0.01])
+            axs[2][1].grid()
 
 
-        # plt.xlabel('Frequency bin')
-        # plt.ylabel('Magnitude |.|^2')
-        f.suptitle('STI output for beam ' + str(beam_to_plot))
+        # # plt.xlabel('Frequency bin')
+        # # plt.ylabel('Magnitude |.|^2')
+        # f.suptitle('STI output for beam ' + str(beam_to_plot))
         # plt.grid(True)
-        plt.savefig("pfb.jpg")
+        plt.savefig("pfbsti_skdjh.png")
 
 except Exception as e:
      print(e)
 
-# bf_len = ntime * nbin * nbeam
+# 
 
 # # RAW BF MODE
-
-# try:
-#     # print("Testing raw mode")
-#     input_len = bf_len * nblocks * 2
+plt.clf()
+try:
+    print("Testing raw mode")
+    input_len = bf_len * nblocks * 2
     
-#     with open("beamformed_6_v2.bin", mode='rb') as file: # b is important -> binary
-#         input_data_f2 = np.array(struct.unpack("f" * (input_len), file.read(4*input_len)))
-#         input_data_c2 = np.array([ complex(input_data_f2[i], input_data_f2[i+1]) for i in range(0,input_len,2) ] ,dtype=complex)
-#     # print(input_data_c2.shape)
-#     stride = nbeam
+    with open("beamformed_1_v2.bin", mode='rb') as file: # b is important -> binary
+        input_data_f2 = np.array(struct.unpack("f" * (input_len), file.read(4*input_len)))
+        input_data_c2 = np.array([ complex(input_data_f2[i], input_data_f2[i+1]) for i in range(0,input_len,2) ] ,dtype=complex)
+    # print(input_data_c2.shape)
+    stride = nbeam
 
-#     print(np.where(input_data_c2 == complex(0,0))[0])
+    print("ZEROS:",np.where(input_data_c2 == complex(0,0))[0])
+    
     
-#     if plot_on:
-#         for c in range(0,nbin):
-#             offset = beam_to_plot + ntime*nbeam*c + block_to_plot*bf_len
-#             indexes = np.arange(offset,offset+ntime*stride,stride,dtype=int)
-
-#             time_in_data_v2 = input_data_c2[indexes]
-#             # print(input_data_c2.shape, time_in_data_v2.shape, offset, stride, bf_len*(block_to_plot+1))
-#             f, input_fft2 = signal.welch(time_in_data_v2, fs, detrend=False, scaling='spectrum', nperseg=512, return_onesided=False)
-#             print(c, max(input_fft2))
-#             plt.plot((fs/2 + fs*c + f)/1e6, 10*np.log10((abs(input_fft2))),'g')
+    if plot_on:
+        for c in range(0,nbin):
+            offset = beam_to_plot + ntime*nbeam*c + block_to_plot*bf_len
+            indexes = np.arange(offset,offset+ntime*stride,stride,dtype=int)
+
+            time_in_data_v2 = input_data_c2[indexes]
+            # print(time_in_data_v2[0:255*9:255])
+            # print(input_data_c2.shape, time_in_data_v2.shape, offset, stride, bf_len*(block_to_plot+1))
+            f, input_fft2 = signal.welch(time_in_data_v2, fs, detrend=False, scaling='spectrum', nperseg=512, return_onesided=False)
+            print(c, max(input_fft2))
+            plt.plot((fs/2 + fs*c + f)/1e6, 10*np.log10((abs(input_fft2))),'g')
             
-#             # print(c, max(abs(time_in_data_v2-time_in_data_v1)))
-#         plt.xlabel('Frequency (MHz)')
-#         plt.ylabel('Power (dB)')
-#         plt.title('Input data FFT')
-
-#         plt.grid(True)
-#         plt.savefig("bf_6" + ".png")
-#         plt.clf()
-# except Exception as e:
-#     print("not testing raw")
-#     print(e)
+            # print(c, max(abs(time_in_data_v2-time_in_data_v1)))
+        plt.xlabel('Frequency (MHz)')
+        plt.ylabel('Power (dB)')
+        plt.title('Input data FFT')
+
+        plt.grid(True)
+        plt.savefig("bf_6" + ".png")
+        plt.clf()
+except Exception as e:
+    print("not testing raw")
+    print(e)
 
 
 

pretty_plot.py

+# TODO:
+# sti and psd are structured the way I think I want everything
+# meaning I didn't change the structure
+# also, you don't need 2 for loops per mode...
+# how to plot sti and psd?
+
+import sys
+import struct
+import matplotlib.pyplot as plt 
+import matplotlib.ticker as ticker
+from scipy import signal, fft
+import numpy as np
+import os
+
+nblocks = 2 #number of blocks saved to each bin file
+plot_on = True
+
+ntime = 5376
+nbeam = 80
+nbin = 25
+fs = 500e6 
+beam_to_plot = 0
+block_to_plot = 1
+bf_len = ntime * nbin * nbeam
+nfft = 256
+
+f, axs = plt.subplots(3,3) #,  constrained_layout = True
+
+try:
+    out_len_streams = nfft * nbin * nbeam//2 * 4 #4 floats / stivec
+    read_out_len = nblocks * out_len_streams
+
+    with open("out_6_v2_c4.bin", mode='rb') as file: # b is important -> binary
+        output_data_f1 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+        output_data_0 = output_data_f1[0:read_out_len:4]
+        output_data_1 = output_data_f1[1:read_out_len:4]
+        output_data_c2 = np.array([ complex(output_data_f1[i], output_data_f1[i+1]) for i in range(2,read_out_len,4) ] )
+
+    if plot_on: 
+        stride = nbeam/2*nbin
+        for c in range(0,1):
+            offset = beam_to_plot + c*nbeam//2 + block_to_plot*len(output_data_0)//2
+            indexes = np.arange(offset,offset+nfft*stride,stride,dtype=int)
+            
+            time_in_data = output_data_0[indexes]
+            axs[2][0].plot(np.arange(0,nfft)+c*nfft, (time_in_data))#, 'C'+str(c % 2)+'o')
+            axs[2][0].set_title("XX* (4 blks)")
+            axs[2][0].set_ylabel("Magnitude")
+            axs[2][0].set_xlabel("Frequency bin")
+            axs[2][0].grid()
+
+            time_in_data = output_data_1[indexes]
+            axs[2][1].plot(np.arange(0,nfft)+c*nfft, (time_in_data))#, 'C'+str(c % 2)+'o')
+            axs[2][1].set_title("YY* (4 blks)")
+            axs[2][1].set_ylabel("Magnitude")
+            axs[2][1].set_xlabel("Frequency bin")
+            axs[2][1].grid()
+
+            
+            time_in_data = output_data_c2[indexes]
+            axs[2][2].plot(np.arange(0,nfft)+c*nfft, (np.abs(time_in_data)))#, 'C'+str(c % 2)+'o')
+            axs[2][2].set_title("XY*+jYX* (4 blks)")
+            axs[2][2].set_ylabel("Magnitude")
+            axs[2][2].set_xlabel("Frequency bin")
+            axs[2][2].grid()
+
+except Exception as e:
+    print(e)
+    plt.clf()
+    f, axs = plt.subplots(2,3)
+    
+
+# K=0.25
+# ntime_plot=ntime
+# x_re = np.round(16*np.random.normal(0,1,ntime_plot) + np.cos(2*np.pi*K*np.arange(0,ntime_plot)))
+# x_im = np.round(16*np.random.normal(0,1,ntime_plot) + np.cos(2*np.pi*K*np.arange(0,ntime_plot)))
+# x = np.array([ complex(x_re[i], x_im[i+1]) for i in range(0,ntime_plot-1) ] ,dtype=complex)
+
+
+# axs[0][0].plot( 10*np.log10((fft.fft(abs(x)))) )
+# axs[0][0].set_xlabel('FFT index')
+# axs[0][0].set_ylabel('Magnitude')
+# axs[0][0].set_title('Beamformer input (FFT)')
+# axs[0][0].grid()
+
+# axs[0][0].plot( np.round(16*np.random.normal(0,1,ntime_plot) + np.cos(2*np.pi*K*np.arange(0,ntime_plot))) )
+# axs[0][0].set_xlabel('Time index')
+# axs[0][0].set_ylabel('Magnitude')
+# axs[0][0].set_title('Beamformer input (real)')
+# axs[0][0].grid()
+
+# axs[0][1].plot( np.round(16*np.random.normal(0,1,ntime_plot) + np.cos(2*np.pi*K*np.arange(0,ntime_plot))) )
+# axs[0][1].set_xlabel('Time index')
+# axs[0][1].set_ylabel('Magnitude')
+# axs[0][1].set_title('Beamformer input (imag)')
+# axs[0][1].grid()
+
+times = [ 0, 10, 20] #[0,6,8]
+try:
+    read_out_len = bf_len*2*nblocks
+    
+    with open("pfb_6_v2.bin",mode='rb') as file:
+        output_data_f2 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+        input_data_c2 = np.array([ complex(output_data_f2[i], output_data_f2[i+1]) for i in range(0,read_out_len,2) ] ,dtype=complex)
+        stride = nbin*nbeam
+        print(output_data_f2[0],input_data_c2[0])
+        print(np.where(input_data_c2 == complex(0,0))[0])
+        if plot_on:
+            
+            for t in range(0,len(times)):
+
+                for c in range(0,1):
+                    offset = beam_to_plot + nbeam*c + times[t]*nbeam*nbin*nfft + block_to_plot*bf_len
+                    indexes = np.arange(offset,offset+nfft*stride,stride,dtype=int)
+
+                    time_in_data_v2 = input_data_c2[indexes]
+                    print(c, max(abs(time_in_data_v2)))
+                    print(abs(time_in_data_v2[0:20]))
+
+                    axs[0][t].plot((abs(time_in_data_v2)))
+                    axs[0][t].set_xlabel('Frequency bin')
+                    axs[0][t].set_ylabel('Magnitude')
+                    axs[0][t].set_title('PFB output t(' + str(times[t]) + ')')
+                    axs[0][t].grid()
+                
+                
+except Exception as e:
+    print("pfb plotting went wrong")
+    print(e)
+
+
+f.tight_layout()
+
+
+try:
+    
+    out_len_streams = nfft * nbin * nbeam//2 * 4 #4 floats / stivec
+    read_out_len = nblocks * out_len_streams
+
+    stride = nbeam/2*nbin
+    with open("out_6_v2_c1.bin", mode='rb') as file: # b is important -> binary
+        output_data_f1 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+        print("sdfklhsdf")
+        output_data_0 = output_data_f1[0:read_out_len:4]
+        print("a;sldfh", len(output_data_0))
+        output_data_1 = output_data_f1[1:read_out_len:4]
+        print("sdjkgf", len(output_data_1))
+        output_data_c2 = np.array([ complex(output_data_f1[i], output_data_f1[i+1]) for i in range(2,read_out_len,4) ] )
+        print("quipr", len(output_data_c2))
+
+    if plot_on: 
+       
+        # f.set_size_inches(6, 6)
+        for c in range(0,1):
+            offset = beam_to_plot + c*nbeam//2 + block_to_plot*len(output_data_0)//2
+            indexes = np.arange(offset,offset+nfft*stride,stride,dtype=int)
+            
+            time_in_data = output_data_0[indexes]
+            print(c,time_in_data)
+            axs[1][0].plot(np.arange(0,nfft)+c*nfft, (time_in_data))#, 'C'+str(c % 2)+'o')
+            axs[1][0].set_title("XX* (1 blk)")
+            axs[1][0].set_ylabel("Magnitude")
+            axs[1][0].set_xlabel("Frequency bin")
+            axs[1][0].grid()
+
+            time_in_data = output_data_1[indexes]
+            axs[1][1].plot(np.arange(0,nfft)+c*nfft, (time_in_data))#, 'C'+str(c % 2)+'o')
+            axs[1][1].set_title("YY* (1 blk)")
+            axs[1][1].set_ylabel("Magnitude")
+            axs[1][1].set_xlabel("Frequency bin")
+            axs[1][1].grid()
+
+            time_in_data = output_data_c2[indexes]
+            axs[1][2].plot(np.arange(0,nfft)+c*nfft, (np.abs(time_in_data)))#, 'C'+str(c % 2)+'o')
+            axs[1][2].set_title("XY*+jYX* (1 blk)")
+            axs[1][2].set_ylabel("Magnitude")
+            axs[1][2].set_xlabel("Frequency bin")
+            axs[1][2].grid()
+
+            # for ax in axs:
+            #     ax.set_yscale('log')
+
+        # tick_spacing = 5000
+        # axs.xaxis.set_major_locator(ticker.MultipleLocator(tick_spacing))
+
+        # # plt.xlabel('Frequency bin')
+        # # plt.ylabel('Magnitude |.|^2')
+        # f.suptitle('STI output for beam ' + str(beam_to_plot))
+        # plt.grid(True)
+except Exception as e:
+    print(e)
+
+
+
+f.suptitle('PFB output (beam 0 / channel 0)')
+f.tight_layout()
+print("Saving figure...")
+plt.savefig("pfbsti_exp_noise_" + str(sys.argv[1]) + ".png")
+plt.clf()
+# 
+
+# # RAW BF MODE
+
+# try:
+#     # print("Testing raw mode")
+#     input_len = bf_len * nblocks * 2
+    
+#     with open("beamformed_6_v2.bin", mode='rb') as file: # b is important -> binary
+#         input_data_f2 = np.array(struct.unpack("f" * (input_len), file.read(4*input_len)))
+#         input_data_c2 = np.array([ complex(input_data_f2[i], input_data_f2[i+1]) for i in range(0,input_len,2) ] ,dtype=complex)
+#     # print(input_data_c2.shape)
+#     stride = nbeam
+
+#     print(np.where(input_data_c2 == complex(0,0))[0])
+    
+#     if plot_on:
+#         for c in range(0,3):
+#             offset = beam_to_plot + ntime*nbeam*c + block_to_plot*bf_len
+#             indexes = np.arange(offset,offset+ntime*stride,stride,dtype=int)
+
+#             time_in_data_v2 = input_data_c2[indexes]
+#             # print(input_data_c2.shape, time_in_data_v2.shape, offset, stride, bf_len*(block_to_plot+1))
+#             f, input_fft2 = signal.welch(time_in_data_v2, fs, detrend=False, scaling='spectrum', nperseg=512, return_onesided=False)
+#             print(c, max(input_fft2))
+#             # axs[1][c].plot((fs/2 + fs*c + f)/1e6, 10*np.log10((abs(input_fft2))),'g')
+            
+#             # print(c, max(abs(time_in_data_v2-time_in_data_v1)))
+#         plt.xlabel('Frequency (MHz)')
+#         plt.ylabel('Power (dB)')
+#         plt.title('Input data FFT')
+
+#         plt.grid(True)
+#         # plt.savefig("bf_6" + ".png")
+#         # plt.clf()
+# except Exception as e:
+#     print("not testing raw")
+#     print(e)
+
+
+
+# try:
+#     # times = [0, 1, 2]
+#     print("transposed data")
+#     out_len_streams = ntime * nbin * nbeam *2
+#     read_out_len = nblocks * out_len_streams
+#     stride = nbeam * nbin
+#     with open("T_6_v2.bin", mode='rb') as file: # b is important -> binary
+#         output_data_f2 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+#         output_data_c2 = np.array([ complex(output_data_f2[i], output_data_f2[i+1]) for i in range(out_len_streams*1,out_len_streams*2,2) ], dtype=complex)
+
+#         if plot_on: 
+#             for c in range(0,1):
+#                 offset = beam_to_plot + c*nbeam
+#                 indexes = np.arange(offset,offset+ntime*stride,stride,dtype=int)
+#                 time_in_data = output_data_c2[indexes]
+               
+#                 # f, input_fft = signal.welch(time_in_data, fs, detrend=False, scaling='spectrum', nperseg=512, return_onesided=False)
+#                 # print(c, 10*np.log10(abs(max(input_fft))))
+#                 # plt.plot((fs/2 + fs*c + f)/1e6, 10*np.log10((abs(input_fft))))
+#                 # print(times)
+#                 print(len(time_in_data),np.where(time_in_data != complex(0,0)))
+
+#                 for t in range(0,len(times)):
+#                     timeDD = np.zeros((nfft,),dtype=complex)
+#                     for i in range(max(0,times[t]-8+1),times[t]+1):
+#                         print(i)
+#                         # print(t, i,time_in_data[nfft*times[t]:nfft*(times[t]+1)])
+#                         timeDD = timeDD + time_in_data[nfft*times[t]:nfft*(times[t]+1)]
+#                     axs[3][t].plot((abs(fft.fft(timeDD))))#np.ones((256,))) #10*np.log10
+#                     axs[3][t].set_xlabel('Frequency bin')
+#                     axs[3][t].set_ylabel('Magnitude')
+#                     axs[3][t].set_title('T6 output t(' + str(times[t]) + ')')
+#                     axs[3][t].grid()
+
+#                     # axs[2][t].plot(abs(timeDD))#np.ones((256,))) #10*np.log10
+#                     # axs[2][t].set_xlabel('time')
+#                     # axs[2][t].set_ylabel('Magnitude')
+#                     # axs[2][t].set_title('T output t(' + str(times[t]) + ')')
+#                     # axs[2][t].grid()
+
+
+#     # print(len(output_data_c2),np.where(output_data_c2 == complex(0,0)))
+#     # plt.plot(np.where(output_data_c2 == complex(0,0))[0])
+#     # plt.savefig("zeros.jpg")
+#     # plt.clf()
+# except Exception as e:
+#     print("idk something went wrong")
+#     print(e)
+
+# try:
+#     with open("T_5_v2.bin", mode='rb') as file: # b is important -> binary
+#         output_data_f2_5 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+#         output_data_c2_5 = np.array([ complex(output_data_f2[i], output_data_f2[i+1]) for i in range(out_len_streams*1,out_len_streams*2,2) ], dtype=complex)
+#         print("diff: ", max(abs(output_data_f2 - output_data_f2_5)))
+#         offset = beam_to_plot# + c*nbeam
+#         indexes = np.arange(offset,offset+ntime*stride,stride,dtype=int)
+#         time_in_data = output_data_c2[indexes]
+#         for t in range(0,len(times)):
+#             timeDD = np.zeros((nfft,),dtype=complex)
+#             for i in range(0,times[t]+1):
+#                 # print(t, i,time_in_data[nfft*times[t]:nfft*(times[t]+1)])
+#                 timeDD = timeDD + time_in_data[nfft*times[t]:nfft*(times[t]+1)]
+#             print(10*np.log10(np.max(abs(fft.fft(timeDD)))))
+#             axs[2][t].plot(10*np.log10(abs(fft.fft(timeDD))),'r')#np.ones((256,))) #
+#             axs[2][t].set_xlabel('Frequency bin')
+#             axs[2][t].set_ylabel('Magnitude')
+#             axs[2][t].set_title('T5 output t(' + str(times[t]) + ')')
+#             axs[2][t].grid()
+
+
+# except Exception as e:
+#     print("guess we're not testing grrr mode", e)
+
+# print("HELLO?")
+
+# PFB
+
+# try:
+    
+#     print("Checking PFB output")
+#     out_len_streams = ntime * nbin * nbeam * 2
+
+#     read_out_len = nblocks * out_len_streams
+#     with open("out_6_v2.bin", mode='rb') as file: # b is important -> binary
+#         output_data_f1 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+#         output_data_c1 = np.array([ complex(output_data_f1[i], output_data_f1[i+1]) for i in range(0,out_len_streams,2) ] )
+
+#     # print(np.where(output_data_c1 == complex(0,0))[0])
+#     # i = np.where(output_data_c1 != complex(0,0))[0]
+
+# except Exception as e:
+#     print("PFB DKFHskldjfa;")
+#     print(e)
\ No newline at end of file

pretty_plot_2.py

+# TODO:
+# sti and psd are structured the way I think I want everything
+# meaning I didn't change the structure
+# also, you don't need 2 for loops per mode...
+# how to plot sti and psd?
+
+import sys
+import struct
+import matplotlib.pyplot as plt 
+import matplotlib.ticker as ticker
+from scipy import signal, fft
+import numpy as np
+import os
+
+nblocks = 2 #number of blocks saved to each bin file
+plot_on = True
+
+ntime = 5376
+nbeam = 80
+nbin = 25
+fs = 500e6 
+beam_to_plot = 0
+block_to_plot = 1
+bf_len = ntime * nbin * nbeam
+nfft = 256
+
+f, axs1 = plt.subplots(1,3) #,  constrained_layout = True
+f.suptitle('PFB output (beam 0 / channel 0)')
+
+# K=0.25
+# ntime_plot=ntime
+# x_re = np.round(16*np.random.normal(0,1,ntime_plot) + np.cos(2*np.pi*K*np.arange(0,ntime_plot)))
+# x_im = np.round(16*np.random.normal(0,1,ntime_plot) + np.cos(2*np.pi*K*np.arange(0,ntime_plot)))
+# x = np.array([ complex(x_re[i], x_im[i+1]) for i in range(0,ntime_plot-1) ] ,dtype=complex)
+
+
+# axs[0][0].plot( 10*np.log10((fft.fft(abs(x)))) )
+# axs[0][0].set_xlabel('FFT index')
+# axs[0][0].set_ylabel('Magnitude')
+# axs[0][0].set_title('Beamformer input (FFT)')
+# axs[0][0].grid()
+
+# axs[0][0].plot( np.round(16*np.random.normal(0,1,ntime_plot) + np.cos(2*np.pi*K*np.arange(0,ntime_plot))) )
+# axs[0][0].set_xlabel('Time index')
+# axs[0][0].set_ylabel('Magnitude')
+# axs[0][0].set_title('Beamformer input (real)')
+# axs[0][0].grid()
+
+# axs[0][1].plot( np.round(16*np.random.normal(0,1,ntime_plot) + np.cos(2*np.pi*K*np.arange(0,ntime_plot))) )
+# axs[0][1].set_xlabel('Time index')
+# axs[0][1].set_ylabel('Magnitude')
+# axs[0][1].set_title('Beamformer input (imag)')
+# axs[0][1].grid()
+
+times = [ 0, 10, 20] #[0,6,8]
+try:
+    read_out_len = bf_len*2*nblocks
+    
+    with open("pfb_6_v2.bin",mode='rb') as file:
+        output_data_f2 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+        input_data_c2 = np.array([ complex(output_data_f2[i], output_data_f2[i+1]) for i in range(0,read_out_len,2) ] ,dtype=complex)
+        stride = nbin*nbeam
+        print(output_data_f2[0],input_data_c2[0])
+        print(np.where(input_data_c2 == complex(0,0))[0])
+        if plot_on:
+            
+            for t in range(0,len(times)):
+
+                for c in range(0,1):
+                    offset = beam_to_plot + nbeam*c + times[t]*nbeam*nbin*nfft + block_to_plot*bf_len
+                    indexes = np.arange(offset,offset+nfft*stride,stride,dtype=int)
+
+                    time_in_data_v2 = input_data_c2[indexes]
+                    print(c, max(abs(time_in_data_v2)))
+                    print(abs(time_in_data_v2[0:20]))
+
+                    axs1[t].plot((abs(time_in_data_v2)))
+                    axs1[t].set_xlabel('Frequency bin')
+                    axs1[t].set_ylabel('Magnitude')
+                    axs1[t].set_title('PFB output t(' + str(times[t]) + ')')
+                    axs1[t].grid()
+                
+                
+except Exception as e:
+    print("pfb plotting went wrong")
+    print(e)
+
+f.tight_layout()
+print("Saving figure...")
+plt.savefig("pfb_exp_noise_" + str(sys.argv[1]) + ".png")
+plt.clf()
+
+f, axs = plt.subplots(2,3) #,  constrained_layout = True
+f.suptitle('STI output (beam 0 / channel 0)')
+try:
+    
+    out_len_streams = nfft * nbin * nbeam//2 * 4 #4 floats / stivec
+    read_out_len = nblocks * out_len_streams
+
+    stride = nbeam/2*nbin
+    with open("out_6_v2_c1.bin", mode='rb') as file: # b is important -> binary
+        output_data_f1 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+        print("sdfklhsdf")
+        output_data_0 = output_data_f1[0:read_out_len:4]
+        print("a;sldfh", len(output_data_0))
+        output_data_1 = output_data_f1[1:read_out_len:4]
+        print("sdjkgf", len(output_data_1))
+        output_data_c2 = np.array([ complex(output_data_f1[i], output_data_f1[i+1]) for i in range(2,read_out_len,4) ] )
+        print("quipr", len(output_data_c2))
+
+    if plot_on: 
+       
+        # f.set_size_inches(6, 6)
+        for c in range(0,1):
+            offset = beam_to_plot + c*nbeam//2 + block_to_plot*len(output_data_0)//2
+            indexes = np.arange(offset,offset+nfft*stride,stride,dtype=int)
+            
+            time_in_data = output_data_0[indexes]
+            print(c,time_in_data)
+            axs[0][0].plot(np.arange(0,nfft)+c*nfft, (time_in_data))#, 'C'+str(c % 2)+'o')
+            axs[0][0].set_title("XX* (1 blk)")
+            axs[0][0].set_ylabel("Magnitude")
+            axs[0][0].set_xlabel("Frequency bin")
+            axs[0][0].grid()
+
+            time_in_data = output_data_1[indexes]
+            axs[0][1].plot(np.arange(0,nfft)+c*nfft, (time_in_data))#, 'C'+str(c % 2)+'o')
+            axs[0][1].set_title("YY* (1 blk)")
+            axs[0][1].set_ylabel("Magnitude")
+            axs[0][1].set_xlabel("Frequency bin")
+            axs[0][1].grid()
+
+            time_in_data = output_data_c2[indexes]
+            axs[0][2].plot(np.arange(0,nfft)+c*nfft, (np.abs(time_in_data)))#, 'C'+str(c % 2)+'o')
+            axs[0][2].set_title("XY*+jYX* (1 blk)")
+            axs[0][2].set_ylabel("Magnitude")
+            axs[0][2].set_xlabel("Frequency bin")
+            axs[0][2].grid()
+
+            # for ax in axs:
+            #     ax.set_yscale('log')
+
+        # tick_spacing = 5000
+        # axs.xaxis.set_major_locator(ticker.MultipleLocator(tick_spacing))
+
+        # # plt.xlabel('Frequency bin')
+        # # plt.ylabel('Magnitude |.|^2')
+        # f.suptitle('STI output for beam ' + str(beam_to_plot))
+        # plt.grid(True)
+except Exception as e:
+    print(e)
+
+try:
+    with open("out_6_v2_c4.bin", mode='rb') as file: # b is important -> binary
+        output_data_f1 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+        output_data_0 = output_data_f1[0:read_out_len:4]
+        output_data_1 = output_data_f1[1:read_out_len:4]
+        output_data_c2 = np.array([ complex(output_data_f1[i], output_data_f1[i+1]) for i in range(2,read_out_len,4) ] )
+
+    if plot_on: 
+       
+        for c in range(0,1):
+            print("sdklasjkl")
+            offset = beam_to_plot + c*nbeam//2
+            indexes = np.arange(offset,offset+nfft*stride,stride,dtype=int)
+            print("qeklhzxvjk")
+            time_in_data = output_data_0[indexes]
+            axs[1][0].plot(np.arange(0,nfft)+c*nfft, (time_in_data))#, 'C'+str(c % 2)+'o')
+            axs[1][0].set_title("XX* (4 blks)")
+            axs[1][0].set_ylabel("Magnitude")
+            axs[1][0].set_xlabel("Frequency bin")
+            axs[1][0].grid()
+            print("aeuauidfzvkl")
+            time_in_data = output_data_1[indexes]
+            axs[1][1].plot(np.arange(0,nfft)+c*nfft, (time_in_data))#, 'C'+str(c % 2)+'o')
+            axs[1][1].set_title("YY* (4 blks)")
+            axs[1][1].set_ylabel("Magnitude")
+            axs[1][1].set_xlabel("Frequency bin")
+            axs[1][1].grid()
+
+            time_in_data = output_data_c2[indexes]
+            axs[1][2].plot(np.arange(0,nfft)+c*nfft, (np.abs(time_in_data)))#, 'C'+str(c % 2)+'o')
+            axs[1][2].set_title("XY*+jYX* (4 blks)")
+            axs[1][2].set_ylabel("Magnitude")
+            axs[1][2].set_xlabel("Frequency bin")
+            axs[1][2].grid()
+
+        
+
+except Exception as e:
+    print(e)
+    for ax in axs[1,:]:
+        f.delaxes(ax)
+        
+f.tight_layout()
+print("Saving figure...")
+plt.savefig("sti_exp_noise_" + str(sys.argv[1]) + ".png")
+plt.clf()
+# 
+
+# # RAW BF MODE
+
+# try:
+#     # print("Testing raw mode")
+#     input_len = bf_len * nblocks * 2
+    
+#     with open("beamformed_6_v2.bin", mode='rb') as file: # b is important -> binary
+#         input_data_f2 = np.array(struct.unpack("f" * (input_len), file.read(4*input_len)))
+#         input_data_c2 = np.array([ complex(input_data_f2[i], input_data_f2[i+1]) for i in range(0,input_len,2) ] ,dtype=complex)
+#     # print(input_data_c2.shape)
+#     stride = nbeam
+
+#     print(np.where(input_data_c2 == complex(0,0))[0])
+    
+#     if plot_on:
+#         for c in range(0,3):
+#             offset = beam_to_plot + ntime*nbeam*c + block_to_plot*bf_len
+#             indexes = np.arange(offset,offset+ntime*stride,stride,dtype=int)
+
+#             time_in_data_v2 = input_data_c2[indexes]
+#             # print(input_data_c2.shape, time_in_data_v2.shape, offset, stride, bf_len*(block_to_plot+1))
+#             f, input_fft2 = signal.welch(time_in_data_v2, fs, detrend=False, scaling='spectrum', nperseg=512, return_onesided=False)
+#             print(c, max(input_fft2))
+#             # axs[1][c].plot((fs/2 + fs*c + f)/1e6, 10*np.log10((abs(input_fft2))),'g')
+            
+#             # print(c, max(abs(time_in_data_v2-time_in_data_v1)))
+#         plt.xlabel('Frequency (MHz)')
+#         plt.ylabel('Power (dB)')
+#         plt.title('Input data FFT')
+
+#         plt.grid(True)
+#         # plt.savefig("bf_6" + ".png")
+#         # plt.clf()
+# except Exception as e:
+#     print("not testing raw")
+#     print(e)
+
+
+
+# try:
+#     # times = [0, 1, 2]
+#     print("transposed data")
+#     out_len_streams = ntime * nbin * nbeam *2
+#     read_out_len = nblocks * out_len_streams
+#     stride = nbeam * nbin
+#     with open("T_6_v2.bin", mode='rb') as file: # b is important -> binary
+#         output_data_f2 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+#         output_data_c2 = np.array([ complex(output_data_f2[i], output_data_f2[i+1]) for i in range(out_len_streams*1,out_len_streams*2,2) ], dtype=complex)
+
+#         if plot_on: 
+#             for c in range(0,1):
+#                 offset = beam_to_plot + c*nbeam
+#                 indexes = np.arange(offset,offset+ntime*stride,stride,dtype=int)
+#                 time_in_data = output_data_c2[indexes]
+               
+#                 # f, input_fft = signal.welch(time_in_data, fs, detrend=False, scaling='spectrum', nperseg=512, return_onesided=False)
+#                 # print(c, 10*np.log10(abs(max(input_fft))))
+#                 # plt.plot((fs/2 + fs*c + f)/1e6, 10*np.log10((abs(input_fft))))
+#                 # print(times)
+#                 print(len(time_in_data),np.where(time_in_data != complex(0,0)))
+
+#                 for t in range(0,len(times)):
+#                     timeDD = np.zeros((nfft,),dtype=complex)
+#                     for i in range(max(0,times[t]-8+1),times[t]+1):
+#                         print(i)
+#                         # print(t, i,time_in_data[nfft*times[t]:nfft*(times[t]+1)])
+#                         timeDD = timeDD + time_in_data[nfft*times[t]:nfft*(times[t]+1)]
+#                     axs[3][t].plot((abs(fft.fft(timeDD))))#np.ones((256,))) #10*np.log10
+#                     axs[3][t].set_xlabel('Frequency bin')
+#                     axs[3][t].set_ylabel('Magnitude')
+#                     axs[3][t].set_title('T6 output t(' + str(times[t]) + ')')
+#                     axs[3][t].grid()
+
+#                     # axs[2][t].plot(abs(timeDD))#np.ones((256,))) #10*np.log10
+#                     # axs[2][t].set_xlabel('time')
+#                     # axs[2][t].set_ylabel('Magnitude')
+#                     # axs[2][t].set_title('T output t(' + str(times[t]) + ')')
+#                     # axs[2][t].grid()
+
+
+#     # print(len(output_data_c2),np.where(output_data_c2 == complex(0,0)))
+#     # plt.plot(np.where(output_data_c2 == complex(0,0))[0])
+#     # plt.savefig("zeros.jpg")
+#     # plt.clf()
+# except Exception as e:
+#     print("idk something went wrong")
+#     print(e)
+
+# try:
+#     with open("T_5_v2.bin", mode='rb') as file: # b is important -> binary
+#         output_data_f2_5 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+#         output_data_c2_5 = np.array([ complex(output_data_f2[i], output_data_f2[i+1]) for i in range(out_len_streams*1,out_len_streams*2,2) ], dtype=complex)
+#         print("diff: ", max(abs(output_data_f2 - output_data_f2_5)))
+#         offset = beam_to_plot# + c*nbeam
+#         indexes = np.arange(offset,offset+ntime*stride,stride,dtype=int)
+#         time_in_data = output_data_c2[indexes]
+#         for t in range(0,len(times)):
+#             timeDD = np.zeros((nfft,),dtype=complex)
+#             for i in range(0,times[t]+1):
+#                 # print(t, i,time_in_data[nfft*times[t]:nfft*(times[t]+1)])
+#                 timeDD = timeDD + time_in_data[nfft*times[t]:nfft*(times[t]+1)]
+#             print(10*np.log10(np.max(abs(fft.fft(timeDD)))))
+#             axs[2][t].plot(10*np.log10(abs(fft.fft(timeDD))),'r')#np.ones((256,))) #
+#             axs[2][t].set_xlabel('Frequency bin')
+#             axs[2][t].set_ylabel('Magnitude')
+#             axs[2][t].set_title('T5 output t(' + str(times[t]) + ')')
+#             axs[2][t].grid()
+
+
+# except Exception as e:
+#     print("guess we're not testing grrr mode", e)
+
+# print("HELLO?")
+
+# PFB
+
+# try:
+    
+#     print("Checking PFB output")
+#     out_len_streams = ntime * nbin * nbeam * 2
+
+#     read_out_len = nblocks * out_len_streams
+#     with open("out_6_v2.bin", mode='rb') as file: # b is important -> binary
+#         output_data_f1 = np.array(struct.unpack("f" * read_out_len, file.read(4*read_out_len)))
+#         output_data_c1 = np.array([ complex(output_data_f1[i], output_data_f1[i+1]) for i in range(0,out_len_streams,2) ] )
+
+#     # print(np.where(output_data_c1 == complex(0,0))[0])
+#     # i = np.where(output_data_c1 != complex(0,0))[0]
+
+# except Exception as e:
+#     print("PFB DKFHskldjfa;")
+#     print(e)
\ No newline at end of file

run_test.bash

@@ -8,7 +8,21 @@ meson compile
 meson install
 # scp haymorer@ssh.et.byu.edu:/fsh/haymorer/Documents/misc_scripts/resample/b.txt ./build
 NUM_BLOCKS=5
-./demortbf -m 6 -b4
+# ./demortbf -m 1 -b1
+cd ..
+for T in 150 200
+do
+    rm *.bin 
+    ./build/demortbf -m 6 -b2 -t$T
+    mv out_6_v2.bin out_6_v2_c1.bin
+    # ./build/demortbf -m 6 -b8 -t$T -c4
+    # mv out_6_v2.bin out_6_v2_c4.bin
+    python3 pretty_plot.py ${T} #pfb_plot.py
+done
+
+
+
+# scp pfb_6.png pfbsti.jpg haymorer@ssh.et.byu.edu:/fsh/haymorer/Documents/misc_scripts/fine_sti
 # ./demortbf -m 0 -b15
 # ./democorr -c2 -b5
 # ./democorr -c3 -b5

src/cublas_beamformer.h

@@ -9,12 +9,14 @@
 #define BLOCKS_TO_SAVE 2
 #define BINS_TO_SAVE 25
 
-#define TONE_FREQ 15.0f
+#define TONE_FREQ 100.0f //not used by rtbfdemo
+
 
 //#define BIN_WITH_DATA 2
-//#define UNIT_STEP_TEST 1
-#define IMPULSE_TEST 1
-#define IS_NOISY 1
+// #define UNIT_STEP_TEST 1
+// #define IMPULSE_TEST 1
+#define COMP_EXP 1
+#define IS_NOISY 0
 
 #define ASYNC 0
 #define TEST_LIST 1
@@ -28,7 +30,7 @@
 
 #define nStreams 4
 
-#define NUM_OPS 18
+#define NUM_OPS 19
 
 // configurable parameters
 #ifndef RTBF_NANTENNA
@@ -56,8 +58,8 @@
 #define RTBF_STI_GPU_BLOCK 32
 #endif
 
-#define NFFT (RTBF_NTIME/RTBF_STI_LENGTH) 
-#define pfb_output_idx(b, t, k, f) (b + RTBF_NBEAM*f + RTBF_NBEAM*RTBF_NBIN*k + RTBF_NBEAM*RTBF_NBIN*NFFT*t) 
+// #define NFFT (RTBF_NTIME/RTBF_STI_LENGTH) 
+#define pfb_output_idx(b, t, k, f) (b + RTBF_NBEAM*f + RTBF_NBEAM*RTBF_NBIN*k + RTBF_NBEAM*RTBF_NBIN*PFB_NFFT*t) 
 //(k + b*NFFT + f*RTBF_NBEAM*NFFT  + t*RTBF_NBIN*RTBF_NBEAM*NFFT)
 #define sti_output_idx_fine(p,b,f,k) (p + b*RTBF_NSTOKES_FLOATS + f*RTBF_CROSSPOL_OFFSET*RTBF_NSTOKES_FLOATS + k*RTBF_CROSSPOL_OFFSET*RTBF_NBIN*RTBF_NSTOKES_FLOATS)
 //(p + k*RTBF_NSTOKES_FLOATS + b*NFFT*RTBF_NSTOKES_FLOATS + f*RTBF_CROSSPOL_OFFSET*NFFT*RTBF_NSTOKES_FLOATS)
@@ -210,7 +212,7 @@ typedef struct RTBFContextStruct {
 
   // host memory pointers
   complex16_t* ibuf_array_h;
-  void* obuf_beamformed_h;
+  float* obuf_beamformed_h;
 
   // host memory buffer sizes
   // these must be sized at least as large as a single input/ouput to the rtbf

src/input_generator.h

+#ifndef _INPUT_GEN_H_
+#define _INPUT_GEN_H_
+
+
+float rand_val(int seed) {
+  const long  a =      16807;  // Multiplier
+  const long  m = 2147483647;  // Modulus
+  const long  q =     127773;  // m div a
+  const long  r =       2836;  // m mod a
+  static long x;               // Random int value
+  long        x_div_q;         // x divided by q
+  long        x_mod_q;         // x modulo q
+  long        x_new;           // New x value
+
+  // Set the seed if argument is non-zero and then return zero
+  if (seed > 0) {
+    x = seed;
+    return(0.0);
+  }
+
+  // RNG using integer arithmetic
+  x_div_q = x / q;
+  x_mod_q = x % q;
+  x_new = (a * x_mod_q) - (r * x_div_q);
+  if (x_new > 0)
+    x = x_new;
+  else
+    x = x_new + m;
+
+  // Return a random value between 0.0 and 1.0
+  return((float) x / m);
+}
+
+/*
+ * Function to generate normally distributed random variable using the Box-Muller method
+ *  params:
+ *    double mean, double standard deviation
+ *  returns:
+ *    normally distributed random variable
+ */
+float norm(float mean, float std_dev) {
+  float   u, r, theta;           // Variables for Box-Muller method
+  float   x;                     // Normal(0, 1) rv
+  float   norm_rv;               // The adjusted normal rv
+
+  //Generate u
+  u = 0.0;
+  while (u == 0.0)
+    u = rand_val(0);
+
+  // Compute r
+  r = sqrt(-2.0 * log(u));
+
+  // Generate theta
+  theta = 0.0;
+  while (theta == 0.0)
+    theta = 2.0 * M_PI * rand_val(0);
+
+  // Generate x value
+  x = r * cos(theta);
+  // TODO: extend for complex random numbers (compute sin value)
+
+  // Adjust x value for specified mean and variance
+  norm_rv = (x * std_dev) + mean;
+
+  // Return the normally distributed RV value
+  return(norm_rv);
+}
+
+
+//first restructure input
+//before the first restructure
+void fillT1Sinusoid(RTBFContext* context) {
+  int Nm = context->mcnt_per_data_block;
+  int Nf = context->nfengines;
+  int Nt = context->time_per_mcnt;
+  int Nc = context->freq_per_xid;
+  int Ni = context->input_per_fid;
+  printf("woot woot\n");
+  size_t pfb_in_bytes = sizeof(complex16_t)*compiletime_info.array_len;
+  complex16_t* h_in = (complex16_t*) malloc(pfb_in_bytes);
+  printf("%i\n", pfb_in_bytes);
+  //cudaMemcpy(h_in, ((RTBFInternalContext*)context->internal)->input_d_start, 
+  //  pfb_in_bytes, cudaMemcpyDeviceToHost);
+ 
+  for (int m=0; m<Nm; m++) {         // mcnt
+    for (int f=0; f<Nf; f++) {       // fid
+      for (int t=0; t<Nt; t++) {     // time
+        for (int c=0; c<Nc; c++) {   // frequency
+          for (int i=0; i<Ni; i++) { // elements per fid
+            int midx = m*(Ni*Nc*Nt*Nf);        // number of complex items per mcount
+            int fidx = f*(Ni*Nc*Nt);           // number of complex items per fid
+            int eidx = i + c*Ni + t*(Nc*Ni) + fidx + midx;
+            float T = (float) (Nt*m + t);
+            h_in[eidx].re = roundf(16.0f*(IS_NOISY*norm(0,1) + cosf(2*3.14*0.25*T)));
+            h_in[eidx].im = roundf(16.0f*(IS_NOISY*norm(0,1) + sinf(2*3.14*0.25*T)));
+            if (!m && !f && !t && !c && !i) printf("%i=%d+%dj ", eidx,h_in[eidx].re,h_in[eidx].im);
+          }
+        }
+      }
+    }
+  }
+  cudaMemcpy(((RTBFInternalContext*)context->internal)->d_restruct_in, h_in, 
+    pfb_in_bytes, cudaMemcpyHostToDevice);
+  free(h_in);
+}
+
+//first restructure output
+//after the first restructure, before the beamformer
+void fillbfSinusoid(RTBFContext* context) {
+  int Nm = context->mcnt_per_data_block;
+  int Nf = context->nfengines;
+  int Nt = context->time_per_mcnt;
+  int Nc = context->freq_per_xid;
+  int Ni = context->input_per_fid;
+  printf("woot woot\n");
+  size_t pfb_in_bytes = sizeof(cuComplex)*compiletime_info.array_len;
+  cuComplex* h_in = (cuComplex*) malloc(pfb_in_bytes);
+  printf("%i\n", pfb_in_bytes);
+  //cudaMemcpy(h_in, ((RTBFInternalContext*)context->internal)->input_d_start, 
+  //  pfb_in_bytes, cudaMemcpyDeviceToHost);
+ 
+  for (int c=0; c<Nc; c++) {   // frequency
+    for (int m=0; m<Nm; m++) {         // mcnt
+      for (int t=0; t<Nt; t++) {     // time
+        for (int f=0; f<Nf; f++) {       // fid
+          for (int i=0; i<Ni; i++) { // elements per fid
+            // int midx = m*(Ni*Nc*Nt*Nf);        // number of complex items per mcount
+            // int fidx = f*(Ni*Nc*Nt);           // number of complex items per fid
+            int eidx = i + Ni*f + Nf*Ni*t + Nt*Nf*Ni*m + Nm*Nt*Nf*Ni*c;//i + c*Ni + t*(Nc*Ni) + fidx + midx;
+            float T = (float) (Nt*m + t);
+            h_in[eidx].x = roundf(16.0f*(IS_NOISY*norm(0,1) + cosf(2*3.14*0.25*T)));
+            h_in[eidx].y = roundf(16.0f*(IS_NOISY*norm(0,1) + sinf(2*3.14*0.25*T)));
+          }
+        }
+      }
+    }
+  }
+  cudaMemcpy(((RTBFInternalContext*)context->internal)->d_array_data, h_in, 
+    pfb_in_bytes, cudaMemcpyHostToDevice);
+  free(h_in);
+}
+
+//beamformer output
+//after the beamformer, before the second restructure
+void fillTSinusoid(RTBFContext* context) {
+  printf("woot woot\n");
+  size_t pfb_in_bytes = sizeof(cuComplex)*compiletime_info.beamformed_len;
+  cuComplex* h_in = (cuComplex*) malloc(pfb_in_bytes);
+  //cudaMemcpy(h_in, ((RTBFInternalContext*)context->internal)->input_d_start, 
+  //  pfb_in_bytes, cudaMemcpyDeviceToHost);
+  int idx = 0;
+  for (int c = 0; c < compiletime_info.nbin; c++) {
+    for (int t = 0; t < compiletime_info.ntime; t++) {
+      for (int b = 0; b < compiletime_info.nbeam; b++) {
+        // printf("%i \n", idx);
+        h_in[idx].x = roundf(9.0f*(IS_NOISY*norm(0,1) + cosf(2*3.14*0.25*t)));
+        h_in[idx].y = roundf(9.0f*(IS_NOISY*norm(0,1) + sinf(2*3.14*0.25*t)));
+
+        // h_in[idx].x = 16.0 * roundf(cosf(2*3.14*0.25*t));
+        // h_in[idx].y = 16.0 * roundf(sinf(2*3.14*0.25*t));
+
+        // if (t == 0) {
+        //   h_in[idx].x = 100.0;
+        //   h_in[idx].y = 0.0;
+        // }
+        // else {
+        //   h_in[idx].x = 0.0;
+        //   h_in[idx].y = 0.0;
+        // }
+
+        // h_in[idx].x = 1.0;
+        // h_in[idx].y = 0.0;
+
+        idx++;
+      }
+    }
+  }
+  cudaMemcpy(((RTBFInternalContext*)context->internal)->d_beamformed, h_in, 
+    pfb_in_bytes, cudaMemcpyHostToDevice);
+  free(h_in);
+}
+
+//second restructure output
+//post second restructure, pre PFB data
+void fillPFBSinusoid(RTBFContext* context) {
+  size_t pfb_in_bytes = sizeof(cuComplex)*compiletime_info.beamformed_len;
+  cuComplex* h_in = (cuComplex*) malloc(pfb_in_bytes);
+  //cudaMemcpy(h_in, ((RTBFInternalContext*)context->internal)->input_d_start, 
+  //  pfb_in_bytes, cudaMemcpyDeviceToHost);
+  int idx = 0;
+  for (int t = 0; t < compiletime_info.ntime; t++) {
+    for (int c = 0; c < compiletime_info.nbin; c++) {
+      for (int b = 0; b < compiletime_info.nbeam; b++) {
+        // printf("%i \n", idx);
+        // h_in[idx].x = roundf(9.0f*IS_NOISY*norm(0,1) + cosf(2*3.14*0.25*t));
+        // h_in[idx].y = roundf(9.0f*IS_NOISY*norm(0,1) + sinf(2*3.14*0.25*t));
+
+        h_in[idx].x = 16.0 * roundf(cosf(2*3.14*0.25*t));
+        h_in[idx].y = 16.0 * roundf(sinf(2*3.14*0.25*t));
+
+        // if (t == 0) {
+        //   h_in[idx].x = 100.0;
+        //   h_in[idx].y = 0.0;
+        // }
+        // else {
+        //   h_in[idx].x = 0.0;
+        //   h_in[idx].y = 0.0;
+        // }
+
+        // h_in[idx].x = 1.0;
+        // h_in[idx].y = 0.0;
+
+        idx++;
+      }
+    }
+  }
+  cudaMemcpy(((RTBFInternalContext*)context->internal)->input_d_start, h_in, 
+    pfb_in_bytes, cudaMemcpyHostToDevice);
+  free(h_in);
+}
+
+#endif
\ No newline at end of file

src/linked_list.cu

@@ -35,19 +35,27 @@
 #define GRRR_DOWN_FACTOR 4
 #endif
 
+typedef struct KernelInfoStruct {
+  dim3 gridDim;
+  dim3 blockDim;
+  void *d_in, *d_out;
+  // int extra_params[1]; //arbitrary number of extra parameters
+  void (*kernel_ptr)(void*,void*);
+} KernelInfo;
 // create a node of the linked list
 typedef struct {
   void* next_op_ptr;                    // point to the next node
   int (*op_function)(RTBFContext*);     // function handle to be called
   //char name[5];                         
   int func_id;                          // function id for function handle
+  KernelInfo kernelInfo;
 } node;
 
 // function identifiers for functions that need to be called 
 // these index an array of functions that can be called from the linked list
 // the array of function handles lives in linked_list.h
 enum func_ids {END = -1,MEMCPY_ON,RESTRUC1,BEAMFORM,SAVE_BF,STI,//SAVE_STI,SAVE_PSD,
-  PSD,RESTRUC2,SAVE_T,GRRR,MEMCPY_OFF,CORR,CHECK_DUMP,CLEAR_BUF,PFB,PFB_STI,SAVE_OUT}; //SAVE_RRR,SAVE_CORR,,SAVE_OUT
+  PSD,RESTRUC2,SAVE_T,GRRR,MEMCPY_OFF,CORR,CHECK_DUMP,CLEAR_BUF,PFB,PFB_STI,SAVE_OUT,SAVE_PFB}; //SAVE_RRR,SAVE_CORR,,SAVE_OUT
 
 /*** basic linked list functions ***/
 // add a new node to the linked list after curr_ptr
@@ -59,6 +67,7 @@ node* newNode(node* curr_ptr, int op, int (*func)(RTBFContext*), bool isHead) {
     curr_ptr->func_id = op;
     curr_ptr->op_function = func;
     curr_ptr->next_op_ptr = NULL;
+    
     return curr_ptr;
 }
 
@@ -106,6 +115,7 @@ int removeNode(node* head, int op_remove) {
     return 0;
 }
 
+
 #define NUM_MODES 7
 /** Context manager **/
 typedef struct LinkedListInternalContext {
@@ -119,10 +129,11 @@ typedef struct LinkedListInternalContext {
   node* checkDump; // points to the node that checks context.doDump
   node* d2h_node; //points to the node that performs the memcpy from device to host
 
-  // output buffer for the current mode:
+  // DEVICE output buffer for the current mode:
   float* d_output; //either points to d_beamformed or a new buffer
   int out_length;
   int out_bytes;
+
   float* h_output; //for reading
 
   //output rtbf_info for all modes
@@ -139,12 +150,27 @@ typedef struct LinkedListInternalContext {
 
   size_t totCudaMalloc;
 
+  // KernelInfo* all_kernelInfo;
+
 } LLInternal;
 
+
+int addKernelInfo(node* node, void (*func)(void*,void*),
+                  dim3 gridDim, dim3 blockDim,
+                  void* in, void* out) {
+  if (!node) return 1;
+  node->kernelInfo.gridDim = gridDim;
+  node->kernelInfo.blockDim = blockDim;
+  node->kernelInfo.d_in = in;
+  node->kernelInfo.d_out = out;
+  node->kernelInfo.kernel_ptr = func;
+  return 0;
+}
+
 //this exists solely for debugging purposes
 char* function_names[NUM_OPS] = { "MEMCPY_ON","RESTRUC1","BEAMFORM","SAVE_BF","STI",//"SAVE_STI","SAVE_PSD"
   "PSD","RESTRUC2","SAVE_T","GRRR","MEMCPY_OFF", //"SAVE_GRRR",
-  "CORR","CHECK_DUMP","CLEAR_BUF","PFB","PFB_STI","SAVE_OUT" }; //,"SAVE_OUT""SAVE_CORR",
+  "CORR","CHECK_DUMP","CLEAR_BUF","PFB","PFB_STI","SAVE_OUT", "SAVE_PFB" }; //,"SAVE_OUT""SAVE_CORR",
 
 // print all the function ids in the linked list
 void print_list(node* head) {
@@ -169,10 +195,10 @@ node* create_linked_list(LLInternal* internal, int* op, int op_mode) {
   bool isHead = true;
   node* curr_ptr = head;
   while (op != NULL && *op != END) {
-    curr_ptr = newNode(curr_ptr,*op++,internal->functions[*op],isHead); 
+    curr_ptr = newNode(curr_ptr,*op++,internal->functions[*op],isHead);
     
     isHead = false; //set this false after the first node
-    // printf("Initialized %s(%i)...%p\n",function_names[curr_ptr->func_id],curr_ptr->func_id,curr_ptr->op_function);
+    printf("Initialized %s(%i)...%p\n",function_names[curr_ptr->func_id],curr_ptr->func_id,curr_ptr->op_function);
 
   }
   
@@ -191,6 +217,7 @@ node* create_linked_list(LLInternal* internal, int* op, int op_mode) {
     // }
     // if (op_mode != BEAM_OP_RAW && op_mode != BEAM_OP_COR && op_mode != BEAM_OP_GRRR)
     curr_ptr = insertAfter(head,MEMCPY_OFF,SAVE_OUT,internal->functions[SAVE_OUT]);
+    curr_ptr = insertAfter(head,PFB,SAVE_PFB,internal->functions[SAVE_PFB]);
   }
 
   return head;
@@ -276,11 +303,14 @@ int init_linked_list(RTBFContext* context,
   //TODO: fix the length, it cannot possibly be correct
   internal->all_out_length[BEAM_OP_PFB] = rtbf_info.sti_pfb_len; // pfb_info.transform_len;
   internal->all_out_bytes[BEAM_OP_PFB] = rtbf_info.sti_pfb_len*sizeof(stokesvec_t); //  pfb_info.transform_len*sizeof(cuComplex);  
-  int f_pfb[10] = {MEMCPY_ON,RESTRUC1,BEAMFORM,RESTRUC2,PFB,PFB_STI,CHECK_DUMP,MEMCPY_OFF,CLEAR_BUF,END}; 
+  int f_pfb[11] = {MEMCPY_ON,RESTRUC1,BEAMFORM,RESTRUC2,PFB,PFB_STI,CHECK_DUMP,MEMCPY_OFF,CLEAR_BUF,END};
+                  //{PFB,PFB_STI,CHECK_DUMP,MEMCPY_OFF,CLEAR_BUF,END};
   internal->all_heads[BEAM_OP_PFB] = create_linked_list(internal,f_pfb,BEAM_OP_PFB);
 
   internal->d_output = NULL;
   internal->checkDump = NULL;
+  // internal->all_kernelInfo = (KernelInfo*) malloc(sizeof(KernelInfo)*NUM_OPS);
+  // for (int i = 0; i < NUM_OPS; i++) internal->all_kernelInfo[i] = 0;
 
   context->list_internal = (void*) internal;
 
@@ -361,13 +391,39 @@ void run_linked_list(RTBFContext* context) {
     while (curr_op != NULL) {
         printf("Running %i(%s)...", curr_op->func_id,function_names[curr_op->func_id]);
         func_timer.Start();
-        (*(curr_op->op_function))(context);
+        KernelInfo kernel_info = curr_op->kernelInfo;
+        if (curr_op->op_function) {
+          (*(curr_op->op_function))(context);
+        }
+        else {
+          printf("%ix%ix%i %ix%ix%i %p %p %p\n", kernel_info.gridDim.x,kernel_info.gridDim.y,kernel_info.gridDim.z,
+          kernel_info.blockDim.x,kernel_info.blockDim.y,kernel_info.blockDim.z, kernel_info.kernel_ptr,
+          kernel_info.d_in,kernel_info.d_out);
+          (*(kernel_info.kernel_ptr))<<<kernel_info.gridDim,kernel_info.blockDim>>>(
+            kernel_info.d_in, kernel_info.d_out);
+        }
         func_timer.Stop();
         printf("runtime: %f msecs.\n", func_timer.Elapsed());
         curr_op = (node*) curr_op->next_op_ptr;
     }
 }
 
+int save_host_data(float* h_data_ptr, int bytes, int mode);
+
+int save_host_data_c(complex16_t* h_data_ptr, int bytes, int mode) {
+  char filename[20];
+  sprintf(filename,"in_%i_v2.bin",mode);
+  FILE *f = fopen(filename, "a");
+
+  if (f != NULL) {
+    fwrite(h_data_ptr, 1, bytes, f);
+    printf("Saved %i bytes to %s ([0] %d+%dj)([252000] %d+%dj)\n", bytes, filename, h_data_ptr[0].re,h_data_ptr[0].im,h_data_ptr[252000].re,h_data_ptr[252000].im);
+  }
+  else printf("file did not open.\n");
+  fclose(f);
+  return RTBF_OK;
+}
+
 int memcpy_to_device(RTBFContext* context) {
 
   LLInternal* internal = (LLInternal*) context->list_internal;
@@ -375,31 +431,62 @@ int memcpy_to_device(RTBFContext* context) {
   // reference to host input and output buffers and adjust for current processing offset
   complex16_t *ibuf_hp = context->ibuf_array_h;
   ibuf_hp += context->input_offset;
+  printf("offset:%i\n",context->input_offset);
 
   #if ASYNC 
 
-  int streamSize = internal->rtbf_info.array_len / nStreams;
-  int streamBytes = internal->in_bytes / nStreams;
+    int streamSize = internal->rtbf_info.array_len / nStreams;
+    int streamBytes = internal->in_bytes / nStreams;
 
-  for (int i = 0; i < nStreams; ++i) {
-    int offset = i * streamSize;
-    // copy memory onto device
-    cudaMemcpyAsync(internal->d_input+offset, ibuf_hp+offset, 
-      streamBytes, cudaMemcpyHostToDevice, internal->stream[i]);
-    checkCudaError();
+    for (int i = 0; i < nStreams; ++i) {
+      int offset = i * streamSize;
+      // copy memory onto device
+      cudaMemcpyAsync(internal->d_input+offset, ibuf_hp+offset, 
+        streamBytes, cudaMemcpyHostToDevice, internal->stream[i]);
+      checkCudaError();
 
-  }
+    }
 
   #else
-
+  save_host_data_c(ibuf_hp,internal->in_bytes,internal->op_mode);
   cudaMemcpy(internal->d_input, ibuf_hp, internal->in_bytes, cudaMemcpyHostToDevice);
   checkCudaError();
 
+  // int Nm = context->mcnt_per_data_block;
+  // int Nf = context->nfengines;
+  // int Nt = context->time_per_mcnt;
+  // int Nc = context->freq_per_xid;
+  // int Ni = context->input_per_fid;
+  // printf("woot woot\n");
+  // complex16_t* h_in = (complex16_t*) malloc(internal->in_bytes);
+  // printf("%i\n", internal->in_bytes);
+ 
+  // for (int m=0; m<Nm; m++) {         // mcnt
+  //   for (int f=0; f<Nf; f++) {       // fid
+  //     for (int t=0; t<Nt; t++) {     // time
+  //       for (int c=0; c<Nc; c++) {   // frequency
+  //         for (int i=0; i<Ni; i++) { // elements per fid
+  //           int midx = m*(Ni*Nc*Nt*Nf);        // number of complex items per mcount
+  //           int fidx = f*(Ni*Nc*Nt);           // number of complex items per fid
+  //           int eidx = i + c*Ni + t*(Nc*Ni) + fidx + midx;
+  //           float T = (float) (Nt*m + t);
+  //           h_in[eidx].re = roundf((cosf(2*3.14*128.0/256.0f*T)));
+  //           h_in[eidx].im = roundf((sinf(2*3.14*128.0/256.0f*T)));
+  //           if (eidx == 252000) printf("%i=%d+%dj ", eidx,h_in[eidx].re,h_in[eidx].im); 
+  //         }
+  //       }
+  //     }
+  //   }
+  // }
+  // save_host_data_c(h_in,internal->in_bytes,7);
+  // // cudaMemcpy(internal->d_input, h_in, internal->in_bytes, cudaMemcpyHostToDevice);
+  // free(h_in);
+
   #endif
   return RTBF_OK;
 }
 
-int save_host_data(float* h_data_ptr, int bytes, int mode);
+
 int save_h_output(RTBFContext* context) {
   LLInternal* internal = (LLInternal*) context->list_internal;
   save_host_data(internal->h_output, internal->out_bytes, internal->op_mode);
@@ -411,17 +498,17 @@ int memcpy_to_host(RTBFContext* context)  {
 
   LLInternal* internal = (LLInternal*)context->list_internal;
 
-  void* obuf_hp = context->obuf_beamformed_h;
-  if (context->operational_mode == BEAM_OP_RAW || 
-      context->operational_mode == BEAM_OP_GRRR || 
-      context->operational_mode == BEAM_OP_COR) {
-    obuf_hp = ((cuComplex*) obuf_hp) + context->output_offset;
-  } else if (context->operational_mode == BEAM_OP_STI || 
-             context->operational_mode == BEAM_OP_PFB) {
-    obuf_hp = ((stokesvec_t*) obuf_hp) + context->output_offset;
-  } else { // BEAM_OP_PSD
-    obuf_hp = ((float*) obuf_hp) + context->output_offset;
-  }
+  float* obuf_hp = ((float*)context->obuf_beamformed_h) + context->output_offset;
+  // if (context->operational_mode == BEAM_OP_RAW || 
+  //     context->operational_mode == BEAM_OP_GRRR || 
+  //     context->operational_mode == BEAM_OP_COR) {
+  //   obuf_hp = ((cuComplex*) obuf_hp) + context->output_offset;
+  // } else if (context->operational_mode == BEAM_OP_STI || 
+  //            context->operational_mode == BEAM_OP_PFB) {
+  //   obuf_hp = ((stokesvec_t*) obuf_hp) + context->output_offset;
+  // } else { // BEAM_OP_PSD
+  //   obuf_hp = ((float*) obuf_hp) + context->output_offset;
+  // }
 
   internal->h_output = (float*) obuf_hp;
 
@@ -454,7 +541,7 @@ int memcpy_to_host(RTBFContext* context)  {
   printf("memcpy %i %p\n", internal->out_bytes, internal->d_output);
   cudaMemcpy(obuf_hp, internal->d_output, internal->out_bytes, cudaMemcpyDeviceToHost);
   checkCudaError();
-  printf("[0] %f [1] %f [2] %f [3] %f\n", ((float*)obuf_hp)[0], ((float*)obuf_hp)[1], ((float*)obuf_hp)[2], ((float*)obuf_hp)[3]);
+  printf("[0] %f [1] %f [2] %f [3] %f [4] %f\n", ((float*)obuf_hp)[0], ((float*)obuf_hp)[1], ((float*)obuf_hp)[2], ((float*)obuf_hp)[3], ((float*)obuf_hp)[4]);
   #endif
 
   return RTBF_OK;
@@ -513,9 +600,10 @@ int LLCleanup(LLInternal* internal) {
   if (internal->op_mode != BEAM_OP_COR && internal->d_cor_matrix != NULL) cudaFree(internal->d_cor_matrix);
 
   cudaFree(internal->d_output);
+  // if (internal->all_kernelInfo) free(internal->all_kernelInfo);
 
   checkCudaError();
 
   return RTBF_OK;
 
-}
\ No newline at end of file
+}

src/restructure.h

@@ -21,8 +21,11 @@
 } while (0)
 #endif
 
+//complex16_t*, cuComplex*
 __global__
-void data_restructure(complex16_t* data, cuComplex* data_restruc) {
+void data_restructure(void* in, void* out) {
+  complex16_t* data = (complex16_t*) in;
+  cuComplex* data_restruc = (cuComplex*) out;
   // Transpose data (corner-turn) from network f-engine order to an order that
   // groups all elements together (removes f-engine dimension) with frequency
   // bin as the slowest moving index. Idea being, why not transpose in the GPU,
@@ -91,9 +94,11 @@ void data_restructure_async(complex16_t * data, cuComplex * data_restruc, int st
 
 }
 
-
+//cuComplex*, cuComplex*
 __global__
-void beam_data_restructure_freq_to_time(cuComplex* data_in, cuComplex* data_out) {
+void beam_data_restructure_freq_to_time(void* in, void* out) {
+  cuComplex* data_in = (cuComplex*) in;
+  cuComplex* data_out = (cuComplex*) out;
   int b = threadIdx.x;
 
   int t = blockIdx.x;

src/rtbf_demo.c

@@ -106,13 +106,16 @@ float norm(float mean, float std_dev) {
   return(norm_rv);
 }
 
-void rtbfFillNetDataBlockBoxMullerRand(complex16_t *h, RTBFContext *context) {
+void rtbfFillNetDataBlockBoxMullerRand(complex16_t *h, RTBFContext *context, int tone_freq) {
   int Nm = context->mcnt_per_data_block;
   int Nf = context->nfengines;
   int Nt = context->time_per_mcnt;
   int Nc = context->freq_per_xid;
   int Ni = context->input_per_fid;
   printf("Filling with complex sinusoidal data.\n");
+
+  float K = (context->operational_mode == BEAM_OP_PFB) ? 1.0*tone_freq / PFB_NFFT : 1.0*tone_freq/2048.0f;
+  printf("NFFT:%i TONE:%i K:%f\n",PFB_NFFT,tone_freq,K);
   // fill network ordered block of array data
   for (int m=0; m<Nm; m++) {         // mcnt
     for (int f=0; f<Nf; f++) {       // fid
@@ -123,7 +126,6 @@ void rtbfFillNetDataBlockBoxMullerRand(complex16_t *h, RTBFContext *context) {
             int fidx = f*(Ni*Nc*Nt);           // number of complex items per fid
             int eidx = i + c*Ni + t*(Nc*Ni) + fidx + midx;
             //printf("%u, %u, %u, %u, %u (%u)\n", m, f, t, c, i, eidx);
-            float K = TONE_FREQ/2048.0f;
             float T = (float) (Nt*m + t); // absolute block time (from 0 to ntime=5376)
 
             #if defined(BIN_WITH_DATA)
@@ -147,11 +149,20 @@ void rtbfFillNetDataBlockBoxMullerRand(complex16_t *h, RTBFContext *context) {
             #elif defined(IMPULSE_TEST)
               h[eidx].re = 1.0; 
               h[eidx].im = 0.0; 
+            #elif defined(COMP_EXP)
+              h[eidx].re = roundf(16.0f*(IS_NOISY*norm(0,1)) + cosf(2*3.14*K*T));
+              h[eidx].im = roundf(16.0f*(IS_NOISY*norm(0,1)) + sinf(2*3.14*K*T));
+              // h[eidx].re = roundf(16.0f*(norm(0,1)) + cosf(2*M_PI*K*T));
+              // h[eidx].im = roundf(16.0f*(norm(0,1)) + sinf(2*M_PI*K*T));
+              // h[eidx].re = roundf((norm(0,1) +cosf(2*3.14*100.0f/2048.0f*T)));
+              // h[eidx].im = roundf((norm(0,1) +sinf(2*3.14*100.0f/2048.0f*T)));
+              // h[eidx].re = roundf(16.0f*(0.2 + cosf(2*3.14*0.25*T)));
+              // h[eidx].im = roundf(16.0f*(0.2 + sinf(2*3.14*0.25*T)));
             #else
-              h[eidx].re = roundf(16.0f*(IS_NOISY*norm(0,1) + cosf(2*M_PI*K*T)));
-              h[eidx].im = roundf(16.0f*(IS_NOISY*norm(0,1) + sinf(2*M_PI*K*T)));
+              h[eidx].re = roundf(16.0f*(IS_NOISY*norm(0,1)));
+              h[eidx].im = roundf(16.0f*(IS_NOISY*norm(0,1)));
             #endif
-            if (!m && !f && !t && !c && !i) printf("%i=%d+%dj ", eidx,h[eidx].re,h[eidx].im);
+            if (eidx == 252000) printf("%i=%d+%dj ", eidx,h[eidx].re,h[eidx].im); //!m && !f && !t && !c && !i
           }
         }
       }
@@ -326,38 +337,37 @@ int rtbfLoadWeightsULABeamResponse(RTBFContext *context, RTBFInfo *rtbfinfo) {
   return rtbf_status;
 }
 
-void* getObuf(RTBFContext* context, RTBFInfo* rtbf_info, PFBInfo* pfb_info, int blocks) {
+void getObuf(RTBFContext* context, RTBFInfo* rtbf_info, PFBInfo* pfb_info, int blocks) {
 
   void *rtbf_out_ptr = NULL;
   size_t output_len;
   switch (context->operational_mode) {
     case BEAM_OP_RAW:
-      rtbf_out_ptr = (float*) malloc(blocks*rtbf_info->beamformed_len*2*sizeof(float));
+      // rtbf_out_ptr = (float*) malloc(blocks*rtbf_info->beamformed_len*2*sizeof(float));
       output_len = blocks*rtbf_info->beamformed_len*2*sizeof(float);
       break;
     case BEAM_OP_STI:
-      rtbf_out_ptr = (stokesvec_t*) malloc(blocks*rtbf_info->sti_stokes_len*sizeof(stokesvec_t));
+      // rtbf_out_ptr = (stokesvec_t*) malloc(blocks*rtbf_info->sti_stokes_len*sizeof(stokesvec_t));
       output_len = blocks*rtbf_info->sti_stokes_len*sizeof(stokesvec_t);
       break;
     case BEAM_OP_PSD:
-      rtbf_out_ptr = (float*) malloc(blocks*rtbf_info->beam_power_len*sizeof(float));
+      // rtbf_out_ptr = (float*) malloc(blocks*rtbf_info->beam_power_len*sizeof(float));
       output_len = blocks*rtbf_info->beam_power_len*sizeof(float);
       break;
     case BEAM_OP_GRRR:
-      rtbf_out_ptr = (float*) malloc(blocks*rtbf_info->resampled_len*2*sizeof(float));
+      // rtbf_out_ptr = (float*) malloc(blocks*rtbf_info->resampled_len*2*sizeof(float));
       output_len = blocks*rtbf_info->resampled_len*2*sizeof(float);
       break;
     case BEAM_OP_PFB:
-      rtbf_out_ptr = (stokesvec_t*) malloc(blocks*rtbf_info->sti_pfb_len*sizeof(stokesvec_t));
+      // rtbf_out_ptr = (stokesvec_t*) malloc(blocks*rtbf_info->sti_pfb_len*sizeof(stokesvec_t));
       output_len = blocks*rtbf_info->sti_pfb_len*sizeof(stokesvec_t);
       break;  
     default:
       printf("error initializing rtbf, unrecognized operational mode\n");
       return NULL;
   }
-  context->obuf_beamformed_h = rtbf_out_ptr; 
+  // context->obuf_beamformed_h = rtbf_out_ptr; 
   context->obuf_beamformed_len = output_len;
-  return rtbf_out_ptr;
 
 }
 
@@ -367,12 +377,13 @@ int main(int argc, char* argv[]) {
   int mode = 0;
   int blocks = 1;
   int dumpOn = 1;
+  int test_tone = 0;
   char weight_filename[PATHLEN];
   int use_weightfile = 0;
   weight_filename[0] = '\0';
   RTBFContext context;
 
-  while ((opt = getopt(argc, argv, "m:b:c:w:vh")) != -1) {
+  while ((opt = getopt(argc, argv, "m:b:c:w:t:vh")) != -1) {
     switch (opt) {
       case 'm':
         // Set operational mode
@@ -409,6 +420,13 @@ int main(int argc, char* argv[]) {
         // be verbose in checking output
         verbose = 1;
         break;
+      case 't':
+          test_tone = strtoul(optarg, NULL, 0);
+          if(test_tone < 0) {
+            fprintf(stderr, "mode must be zero, or greater\n");
+            return 1;
+          }
+        break;
       default:
         fprintf(stderr,
             "Usage: %s [options]\n"
@@ -484,7 +502,11 @@ int main(int argc, char* argv[]) {
   #if TEST_MODES
   void* rtbf_out_ptr = getObuf(&context,&rtbf_info,&pfb_info, blocks < 6 ? blocks : 5 );
   #else
-  void* rtbf_out_ptr = getObuf(&context,&rtbf_info,&pfb_info, blocks);
+  context.obuf_beamformed_len = blocks*rtbf_info.beamformed_len*2*sizeof(float);
+  float* rtbf_out_ptr = (float*) malloc(context.obuf_beamformed_len);
+  context.obuf_beamformed_h = rtbf_out_ptr; 
+  
+  //getObuf(&context,&rtbf_info,&pfb_info, blocks);
   #endif
   // void* rtbf_out_ptr = NULL;
 
@@ -573,8 +595,9 @@ int main(int argc, char* argv[]) {
 
       updateOpMode(&context,mode);
       printf("Changing modes in rtbf_demo! %i\n",mode);
-      free(rtbf_out_ptr);
-      rtbf_out_ptr = getObuf(&context,&rtbf_info,&pfb_info,5);
+      //free(rtbf_out_ptr);
+      //rtbf_out_ptr = 
+      getObuf(&context,&rtbf_info,&pfb_info,5);
       registerBuffers(&context,0);
       rand_val(mode); //reset the seed
     }
@@ -590,34 +613,36 @@ int main(int argc, char* argv[]) {
     b_idx %= 5;
     #endif
 
-    if (mode == BEAM_OP_RAW) {
-      context.output_offset = b_idx*rtbf_info.beamformed_len;
-      // the 2 has to be here and not on setting context.output_offset because
-      // the internal library is doing pointer addition on a cuComplex* but here
-      // is done on a float*
-      update_b_out_ptr = ((float*) rtbf_out_ptr) + 2*context.output_offset;
-    } else if (mode == BEAM_OP_STI) {
-      context.output_offset = b_idx*rtbf_info.sti_stokes_len;
-      update_b_out_ptr = ((stokesvec_t*) rtbf_out_ptr) + context.output_offset;
-    }
-    else if (mode == BEAM_OP_GRRR) {
-      context.output_offset = b_idx*rtbf_info.resampled_len;
-      // the 2 has to be here and not on setting context.output_offset because
-      // the internal library is doing pointer addition on a cuComplex* but here
-      // is done on a float*
-      update_b_out_ptr = ((float*) rtbf_out_ptr) + 2*context.output_offset;
-    }
-    else if (mode == BEAM_OP_PFB) {
-      context.output_offset = b_idx*rtbf_info.sti_pfb_len;
-      // the 2 has to be here and not on setting context.output_offset because
-      // the internal library is doing pointer addition on a cuComplex* but here
-      // is done on a float*
-      update_b_out_ptr = ((stokesvec_t*) rtbf_out_ptr) + context.output_offset;
-    }  
-    else { // BEAM_OP_PSD
-      context.output_offset = b_idx*rtbf_info.beam_power_len;
-      update_b_out_ptr = ((float*) rtbf_out_ptr) + context.output_offset;
-    }
+    context.output_offset = b_idx*rtbf_info.beamformed_len;
+    update_b_out_ptr = ((float*) rtbf_out_ptr) + 2*context.output_offset;
+    // if (mode == BEAM_OP_RAW) {
+    //   context.output_offset = b_idx*rtbf_info.beamformed_len;
+    //   // the 2 has to be here and not on setting context.output_offset because
+    //   // the internal library is doing pointer addition on a cuComplex* but here
+    //   // is done on a float*
+    //   update_b_out_ptr = ((float*) rtbf_out_ptr) + 2*context.output_offset;
+    // } else if (mode == BEAM_OP_STI) {
+    //   context.output_offset = b_idx*rtbf_info.sti_stokes_len;
+    //   update_b_out_ptr = ((stokesvec_t*) rtbf_out_ptr) + context.output_offset;
+    // }
+    // else if (mode == BEAM_OP_GRRR) {
+    //   context.output_offset = b_idx*rtbf_info.resampled_len;
+    //   // the 2 has to be here and not on setting context.output_offset because
+    //   // the internal library is doing pointer addition on a cuComplex* but here
+    //   // is done on a float*
+    //   update_b_out_ptr = ((float*) rtbf_out_ptr) + 2*context.output_offset;
+    // }
+    // else if (mode == BEAM_OP_PFB) {
+    //   context.output_offset = b_idx*rtbf_info.sti_pfb_len;
+    //   // the 2 has to be here and not on setting context.output_offset because
+    //   // the internal library is doing pointer addition on a cuComplex* but here
+    //   // is done on a float*
+    //   update_b_out_ptr = ((stokesvec_t*) rtbf_out_ptr) + context.output_offset;
+    // }  
+    // else { // BEAM_OP_PSD
+    //   context.output_offset = b_idx*rtbf_info.beam_power_len;
+    //   update_b_out_ptr = ((float*) rtbf_out_ptr) + context.output_offset;
+    // }
     printf("offset:%i\n",context.output_offset);
 
     // fill data block
@@ -626,7 +651,7 @@ int main(int argc, char* argv[]) {
       rtbfRAWFillDataBlock(update_d_cx_ptr, &context, &rtbf_info);
     } 
     else {
-      rtbfFillNetDataBlockBoxMullerRand(update_d_cx_ptr, &context);
+      rtbfFillNetDataBlockBoxMullerRand(update_d_cx_ptr, &context,test_tone);
     }
 
     doDump = ((b % dumpOn) == 0 || (b==blocks)) ? 1 : 0;

src/save_bin.h

@@ -105,7 +105,7 @@ int save_device_data(float* d_data_ptr, int bytes, char filename[20]) {
     //   idx = ii*floats_to_read;
     //   if (h_data_ptr[idx] < minVal) minVal = h_data_ptr[idx];
     // } 
-    printf("Saved %i bytes to %s ([0] %f)\n", bytes, filename, h_data_ptr[0]);
+    printf("Saved %i bytes to %s ([0] %f [%i] %f)\n", bytes, filename, h_data_ptr[0], RTBF_NBEAM*RTBF_NBIN*256*8, h_data_ptr[RTBF_NBEAM*RTBF_NBIN*256*8]);
   }
   else printf("file did not open.\n");
   fclose(f);

subprojects/pfb.wrap

 [wrap-git]
 
 url = git@gitlab.ras.byu.edu:alpaca/xb-engine/core/pfb
-revision = master
+revision = LL_compatible
 depth = 1

test.bash

+for T in 150 200
+do
+    rm *.bin 
+    ./build/demortbf -m 6 -b2 -t$T
+    mv out_6_v2.bin out_6_v2_c1.bin
+    python3 pretty_plot.py ${T}
+done
\ No newline at end of file