texts: script to check zero crossings to check if microphones phases are aligned

texts:
I have 'n' microphones placed in square fashion and I want to make sure that all the channels must be aligned in time exactly the same if the signal is equidistant from all the 'n' microphones i.e. in the center of the square.

I have written below script to do the difference in zero crossing timings and if the difference is above some precision then print that and bail out.

from scipy.io import wavfile
import numpy as np
import argparse

parser = argparse.ArgumentParser(description='diff the zero crossing of two files')
parser.add_argument('-f1', '--file_name_1', help='provide first file name')
parser.add_argument('-f2', '--file_name_2', help='provide second file name')
parser.add_argument('-p', '--precision', help='precision to compare against', type=float, default=0.0001)

args = parser.parse_args()
print(args)
files = []
files.append(args.file_name_1)
files.append(args.file_name_2)

sampling_rates = []
signals = []
for file in files:
  fs, signal = wavfile.read(file)
  signal = signal / max(abs(signal))                        # scale signal
  sampling_rates.append(fs)
  signals.append(signal)
  assert min(signal) >= -1 and max(signal) <= 1
  print 'fs           ==> ', fs, 'Hz'                       # sampling rate
  print 'len(signal)  ==> ', len(signal), 'samples'

sampsPerMilli = 0
#files should have same sampling rates and length
prev_rate = sampling_rates[0]
for rate in sampling_rates[1:]:
  cur_rate =  rate
  sampsPerMilli = int(rate / 1000)
  if prev_rate != cur_rate:
    print("rates doesn't match %d %d"% (prev_rate, cur_rate))
    exit(0)
  cur_rate = rate

#signal length also should be same
prev_length = len(signals[0])
for signal in signals[1:]:
  cur_length = len(signal)
  if prev_length != cur_length:
    print("length of signals doesn't match for %d %d"% (prev_length, cur_length))
    exit(0)
  cur_length = prev_length

zccs = []
for signal in signals:
  zcc = []
  DC = np.mean(signal)
  newSignal = signal - DC
  for i in range(1, int(len(newSignal))):
    if((newSignal[i] * newSignal[i-1]) < 0):
      #print("crossing at %f seconds"% ((i/sampsPerMilli) * 0.001))
      zcc.append((i/sampsPerMilli) * 0.001)
  zccs.append(zcc)

for a, b in zip(zccs, zccs[1:]):
  if len(a) != len(b):
    print("length doesn't match %d %d"% (len(a), len(b)))
    for c, d in zip(a, b):
      if c-d > args.precision:
        print("precision %f c %f d %f exceeded"% (args.precision, c, d))
        exit(0)

Is there any better approach or can this script be improved?


texts: How to change volume of a PCM 16 bit signed audio?

texts:
I know I can multiply samples, then clip but perceived volume is non-linear for humans.

Can you please help with a formula.


texts: Does $ Ev(\cos(4\pi t)u(t)) $ signal has period?

texts:
According to Oppenheim's signal and systems solution book and other solutions like this

$$
Ev(\cos(4\pi t)u(t)) = 1/2 \times ( \cos(4\pi t)u(t) + \cos(4\pi t)u(-t)) = 1/2(\cos(4\pi t))
$$

$$
\textrm{for }-\infty < t <  \infty
$$

is periodic and its period $1/2$.

But I think it is not equal $ 1/2(\cos 4\pi t)  $ at $ t = 0 $ so it is not periodic because at $ t = 0 $ it is equal to $ 1 $ different from other multiples of period.

I have also curious about the value of $ u(t)+u(-t) $. Is it $0$, $1$, or $2$ ?


texts: How do I remove all white pixels in a modified image which are not present in the original?

texts:
I have the following image

and after thresholding it, I obtained this result:

How can I remove all white pixels which do not belong to the original image?

The desired output should look like this (I manually edited the image):


texts:
I was reading Winder's Analog and Digital Filter Design and the section on Bessel filter.

I was hoping to see a complete derivation of the Bessel filter theory, but Winder's book gives only

\begin{align}
H(s)&=e^{-sT} \tag{$\scriptstyle\text{pure delay}$}\\
H(s)&=e^{-s}=\frac{1}{\sinh(s)+\cosh(s)} \tag{$\scriptstyle\text{normalized pure delay, T=1}$}
\end{align}

and then using series expansion of $\sinh$ and $\cosh$ we're left with

$$H(s)=\frac{a_0}{B_n(s)}$$

where $B_n(s)$ is a Bessel polynomial.

Anyone know where can I get a more exhaustive derivation of the Bessel filters?


texts: How to do simple extrapolation with Fourier transformation?

texts:
I have 1024 sample points, and I would like to do really simple extrapolation using Fourier transformation. First I apply Fast fourier transformation on the data. My first intuition was that I just calculate the inverse fourier transformation on a larger interval. The result of the FFT is a collection of sinusoidal waves, so I expected that it's result is continuous at the end of the original interval.

I calculate the result with this formula:

$$
X[k] = \sum_{n=0}^{N-1} x[n] e^{\frac{i 2 \pi n k}{N}}
$$

The extrapolation at the end of the interval looks like this:

I would like to avoid this "break" in the result.

The python source code is:

import numpy as np
import matplotlib.pyplot as plt

from scipy.fftpack import fft

# generate samples
data = np.zeros(1024)
for k in range(0, 1024):
    data[k] = np.sin(k/20 * np.pi*2)

# apply FFT on samples
waves = fft(data)

# calculate the new data points from the FFT result
res = np.zeros(2048)
for k in range(0, 2048):
    val = 0.0
    for n in range(0, len(waves)):
        # https://dsp.stackexchange.com/a/510/25943
        val += waves[n]*np.exp(1.j * 2*np.pi * n * k / len(waves)) / len(waves)
    res[k] = val.real

#plot the result
plt.plot(res)
plt.show()


texts: Rocketry : Combine two accelerometers to reduce noise?

texts:
I am designing an IMU for an experimental rocket. I'll be using the BN055 9DOF that has sensor fusion - orientation quaternion - and gravity compensation.

My main goal is speed computation, and I was told that using 2 accelerometers in the same time may help reduce the noise (I'll implement a kalman filter to combine data from both accelerometers and integrate).

So my question is : Is that relevant ? And if yes should the 2nd accelerometer be high G (+100G) or low G (16G same as BN055) ?

Thank you in advance for your help


texts: How does one go about choosing a Radar Cross Section for a human when modelling a Pulse micro-Doppler microwave radar?

texts:
So as part of a personal project I'm in the process of simulating and making a fall detector using a 2.7GHz microwave sensor with a range of 8m. I've been using MATLAB's example code to simulate each part of the radar as well as a target and a distance returned which can be found here: https://www.mathworks.com/help/phased/ug/designing-a-basic-monostatic-pulse-radar.html#d123e23746
As well as a radar based on the BISS0001 Chip with data sheet: https://cdn-learn.adafruit.com/assets/assets/000/010/133/original/BISS0001.pdf
This seems to work quite well but one thing I'm not entirely sure on is the design specifications and how I can translate them to my radar. In MATLAB's example they have a max unambiguous range of 5000m, a resolution range of 50m and a radar cross section of 1m^2.
Given I'm not expecting more than 1 target together I figure I don't really need to worry too much about maximum range and resolution range but regardless, since the range of my radar is 8m I've put the max range as 8m and as for the range resolution I've gone for 2m.
I'm not entirely sure that those are correct but given that I'll never have more than 1 target I think it might be OK but please let me know if not.
The main thing is the Radar Cross Section as it's made up of a lot of properties I'm not definite on. Looking online at published data it seems 1m^2 is standard but does this not depend on the strength of my radar instead of just the target itself?


texts: Need to calculate $\frac{E\left[|x(n)|^4\right]}{E\left[x^2(n)\right]}$

texts:
I need to calculate $\frac{E\left[|x(n)|^4\right]}{E\left[x^2(n)\right]}$ in real time, where $x(n)$ is complex.

For $E\left[|x(n)|^4\right]$ I plan to take the absolute, raise it to the power $4$ and do exponential averaging with weights $0.999$ and $0.001$. I have no clue on how to calculate $E\left[x^2(n)\right]$. I know that this gives the total signal power. Will this be equivalent to $E\left[|x(n)|^2\right]$?

Thanks, JK


texts: Interpretation of Convergence plot in terms of mean square error

texts: I am trying to learn various learning algorithms used in training of neural network. When plotting the Mean Square Error in dB with respect to the number of iterations, if the points are around 0 dB then can we say that the curve converges or does the points need to go to minus dB? If the curve goes towards minus dB with the increase in iterations, then what can be inferred? How good or bad is the technique?


texts: Isolate the non blurred part of foucsed image

texts:
I have an image taken by a mobile phone camera with focus on a part of it let say a face or whatever.
I want to know -approximately-  where is the focused area. For example, a bounding box around the focused area.
example :


texts: What is causing my ANC LMS update to diverge?

texts:
I'm trying to implement a simulation of an ANC system with python, using this model here.
My simulation keeps diverging, and I honestly don't know why. I'm using a source for LMS adaptive filter from Mathworks here. When I comment out the LMS update function, the system is stable, and the output is exactly the input. But when I plug in the adaptive filter update, the system starts to diverge. So I thought the problem is with my LMS update. And I implemented the LMS adaptive filter with padasip package. The system is still diverging. Now I honestly don't know how I copied the model from Matlab wrong. Can someone help? Though It's not a syntax error, I pasted my code below for reference. I used all difference equations for every system.
Some values I used in this,
The input I'm using is a white noise file I generated from Matlab, 1dB power and sampling frequency of 16000 Hz.
b for S(z) is [0.5,0.5,-.3,-.3,-.2,-.2,]
b for S_est(z) is [0.466,0.533,-0.257,-0.274,-0.231,-0.175]
b for Main path S(z) is [0.0500,0,0.0200,0,-0.0000,0,-0.1250,0,-0.0500,0,0.0750,0,0.0300]
All I copied from the Simulink simulation of the same system
#ANC Simulation Main

import pyaudio, wave, struct, math
import numpy as np
from matplotlib import pyplot as plt
import padasip as pa

## Variables Setup
mu = 0.1
MAXVALUE = 2**15-1  # Maximum allowed output signal value (because WIDTH = 2)

# Initialization of adaptive weight w
w = np.zeros(13)

# Main Path Filters
order_path = 12
b_path = np.array([0.0500,0,0.0200,0,-0.0000,0,-0.1250,0,-0.0500,0,0.0750,0,0.0300])
x = np.zeros(13)
y_main_path =0

#Adapt filter
order_adapt = 12
y_adapt =np.zeros(13)

filt = pa.filters.FilterLMS(13, mu=mu)

#second path
order_sec_path = 12
b_sec_path = [0.5,0.5,-.3,-.3,-.2,-.2,0,0,0,0,0,0,0]

y_secpath = 0

#est sec path
order_est_sec_path = 12
b_est_sec_path = np.array([0.466,0.533,-0.257,-0.274,-0.231,-0.175,0,0,0,0,0,0,0])

est_sec_out = 0
y_estsecpath = np.zeros(13)

# File names
wavfile = 'matlab_1db_wn.wav'
output_wavfile = 'ANC_Result.wav'

## Read WAV file
wf = wave.open(wavfile,'rb')

CHANNELS        = wf.getnchannels()     # Number of channels
RATE            = wf.getframerate()     # Sampling rate (frames/second)
signal_length   = wf.getnframes()       # Signal length
WIDTH           = wf.getsampwidth()     # Number of bytes per sample

print('The file has %d channel(s).'            % CHANNELS)
print('The frame rate is %d frames/second.'    % RATE)
print('The file has %d frames.'                % signal_length)
print('There are %d bytes per sample.'         % WIDTH)

# Read first BLOCKLEN
binary_data = wf.readframes(1)

## Output WAV file
output_wf = wave.open(output_wavfile, 'w')
output_wf.setframerate(RATE)
output_wf.setsampwidth(WIDTH)
output_wf.setnchannels(CHANNELS)

## Open audio stream
p = pyaudio.PyAudio()
stream = p.open(
   format      = p.get_format_from_width(WIDTH),
   channels    = CHANNELS,
   rate        = RATE,
   input       = False,
   output      = True )

## Main Loop
while len(binary_data) > 0:

   # convert binary data to numbers
   input_block = struct.unpack('h', binary_data)
   input_value = input_block[0]

   x = np.delete(x,-1)
   x = np.insert(x,0,input_value)

   y_main_path = np.dot(b_path,np.transpose(x))

   y_adapt = filt.predict(input_value)

   #y_adapt = np.delete(y_adapt,-1)
   #y_adapt = np.insert(y_adapt,0,adapt_out)

   est_sec_out = np.dot(b_est_sec_path,np.transpose(x))
   y_estsecpath = np.delete(y_estsecpath,-1)
   y_estsecpath = np.insert(y_estsecpath,0,est_sec_out)

   #sec path
   y_secpath = np.dot(b_sec_path,np.transpose(y_adapt))

   output = y_main_path - input_value

   # LMS update
   filt.adapt(y_main_path,input_value)

   output = np.clip(output, -MAXVALUE, MAXVALUE)
   output = output.astype(int)
   binary_date = struct.pack('h', output)
   stream.write(binary_data)
   output_wf.writeframes(binary_data)
   binary_data = wf.readframes(1)

stream.stop_stream()
stream.close()
p.terminate()
# Close wavefiles
wf.close()
output_wf.close()


texts: Turn a decision problem into a linear filter problem

texts:
Let $x[n]$ be a time-series, and two filters:

$A[n] = a_0 x[n] + a_1 x[n-1] + ... a_q x[n-q]$
$B[n] = b_0 x[n] + b_1 x[n-1] + ... b_r x[n-r]$

I think the answer is No, but is there a way to have this problem (sign of 2 filters should be both positive):

$$A[n] > 0 \quad AND \quad  B[n] > 0$$

expressed in terms of the study of the sign of a single linear filter

$$C[n] = c_0 x[n] + c_1 x[n-1] + ... c_s x[n-s]$$

?

If not possible, is there a way to have nearly equivalence (that holds for most values of $n$)?

Application:

A[n] = x[n] - x[n-1]                 # derivative  (or a smoothed version with moving average)
B[n] = x[n] - 2 * x[n-1] + x[n-2]    # acceleration (idem)

I'd like to find when both derivative and acceleration of $x[n]$ are positive by studying the sign of only 1 linear filter.

Or said in another way: is it possible to create a linear filter $C[n] = c_0 x[n] + c_1 x[n-1] + ... c_s x[n-s]$ such that $C[n] > 0$ if and only if derivative and acceleration of $x[n]$ are both positive?


texts: How do real-time convolution plugins process audio so quickly

texts:
Okay, so in Logic Pro I can load up a Space Designer plugin (convolution reverb) with an impulse that's 9.1 seconds long, turn my mic on, and get real-time convolution reverb as the mic records incoming audio.

My question is: How could this be possible? My I/O buffer size is 32 samples. Is my computer capable of performing the convolution of a 401,310 sample impulse with a 32 sample block of audio in under 0.0007 seconds? If so, what kind of black magic do I need to get my own convolution reverb to do this?

I should mention that I am trying to develop real time convolution in C with portaudio. It works great when it's not real-time, but it's way too slow in real-time (the convolution I mentioned above takes 0.8 seconds using frequency-domain convolution which is supposed to be faster for impulses of that size, no?

Anyway any insight would be helpful, as I would really like to understand how this kind of real-time processing is possible.


texts: What is the difference between a Fourier transform and a cosine transform?

texts:
In speech recognition, the front end generally does signal processing to allow feature extraction from the audio stream. A discrete Fourier transform (DFT) is applied twice in this process. The first time is after windowing; after this Mel binning is applied and then another Fourier transform.

I've noticed however, that it is common in speech recognizers (the default front end in CMU Sphinx, for example) to use a discrete cosine transform (DCT) instead of a DFT for the second operation. What is the difference between these two operations? Why would you do DFT the first time and then a DCT the second time?


texts: Is possible to compare to audio samples with cross-correlation

texts: Is it possible to compare two audio samples(meaning find word from one in another) realtime with cross-correlation?


texts:
Is there an accurate way to create an aliased image from the Fourier transform of the original image?

in other words, i have the Fourier coefficients of an image and i want to make down-sample in frequency domain.


texts: Is there a way to do a DFT with low latency and fine frequency resolution at low frequencies?

texts:
I am learning about the DFT in order to design a fundamental pitch detector for the human singing voice. Using the cepstrum technique.

I am concerned about the latency for low frequency detection at <100Hz. Is it always necessary for the DFT algorithm to work with sample length that contains an entire wavelength of the lowest frequency being detected, or can the window be shorter than the wavelength? For example if the singer sings <100 Hz then there is the possibility of the >10ms latency. This latency could get worse with further processing.

I am also concerned about frequency resolution at low frequencies. The standard DFT is inherently linear this causes poorer pitch discrimination at lower frequencies, If frequency index 1 was 50Hz, then the next bins would be 100Hz and 150Hz. The steps are clearly too big for pitch detection. If the index started at a lower frequency to get finer graduations, say 1 Hz steps, then the latency becomes 1 second, if what I said in the previous paragraph is true. Btw 1Hz steps is still pretty poor for pitch detection.

Are there more sophisticated DFT algorithms that address these issues?

I know products like Antares Auto-Tune have pretty fast detection, so this stuff is possible.


texts: Hough Transform in Matlab without the built-in function

texts:
I want to implement Hough transform on image without using inbuilt function. Some papers say that the image is first flipped before applying Hough transform. I do not understand how Matlab is doing it.

I have written the code below, but the H matrix by Matlab and houghMatrix generated by me are not same.

I have used default values of rhoResolution as 0.5 and thetaResolution as 0.5 also. I tried to do it without rotating image but then also not getting same matrices.

RGB = imread('gantrycrane.png');
     I  = rgb2gray(RGB); % convert to intensity
     BW = edge(I,'canny'); % extract edges
     BW = BW(end:-1:1,1:1:end);
     theta = -90:0.5:(90-0.5);
     D = sqrt((size(BW,1) - 1)^2 + (size(BW,2) - 1)^2);
     rho = -ceil(D):0.5:ceil(D);
     houghMatrix = zeros(size(rho,2),size(theta,2));
     for i=1:size(BW,1)
         for j=1:size(BW,2)
             if(BW(i,j)==1)
                  for ii=1:size(theta,2)
                      rho1 = i*cosd(theta(1,ii)) + j*sind(theta(1,ii));
                      rhotemp = rho1 - floor(rho1);
                      if(rhotemp>(0.5+0.25))
                          rho1 = ceil(rho1);
                      elseif(rhotemp>0.25)
                          rho1 = floor(rho1) + 0.5;
                      else
                          rho1 = floor(rho1);
                      end
                      position = (rho1 - rho(1,1))/0.5;
                      houghMatrix(position+1,ii) = houghMatrix(position+1,ii) + 1;
                  end
             end
         end
     end

 [H,T,R] = hough(BW,'RhoResolution',0.5,'Theta',-90:0.5:89.5);

I need to implement this in C code after that, and the functionality should be same as Matlab code. So I need to have a similar matrix like Matlab (if error of +-5 then it is acceptable but this is two totally different matrices).


texts:
QUESTION

Compute the Z-transform of $y[n] = x[n] + 2x[n-1]$. and find the poles and zeros.

I just bombed an interview where I couldn't do this (because I have no grounding in fundamentals and have worked all from cookbooks). I'm not too bummed because I would have been working at something I wasn't prepared for. You get one guess who would ask this kind of question in a first phone interview.

ANSWER (30 minutes post interview):

$x[n-1] = z^{-1}$

$1 + 2z^{-1}$.

Why?

Because the z-transform of an FIR is the coefficients of the impulse response multiplied by the delays.


texts: Conditions on $v \in \mathbb{C}^{n}$ such that the complex to real DFT is sensible

texts:
I'm using FFTW to compute some DFTs, and as a memory space optimization, you can give it a $v \in \mathbb{R}^{n}$ and it will return a $\hat{v} \in \mathbb{C}^{\lfloor n/2 \rfloor + 1}$ via the call fftw_execute_dft_r2c. (To those without context, it's just claiming that since $\overline{\hat{v}_{k}} = \hat{v}_{n-k}$, then you just don't need to care about the top half of the complex vector.)

In addition, you can take $\hat{v}$ back to $\mathbb{R}^{n}$ via fftw_execute_dft_c2r, but of course no one does that, they want to apply some filter to $\hat{v}$ and then call fftw_execute_dft_c2r.

What is the largest class of operations on $\hat{v}$ such that fftw_execute_dft_c2r make sense?

(Clearly, pointwise multiplication by the DFT of another real function is allowed, but it's not clear to me that every bandpass filter produces a sensible output and many other operations seem suspect.)


texts: Phase response of moving average filter -- how to interpret?

texts:
There are many articles on the frequency response of the moving average filter but they all seem to focus on magnitude. However the phase response is intriguing and I find it hard to interpret. The phase appears to "wrap", but it wraps within the [-pi, pi) interval rather than at its edges. Example below:

A phase unwrapping algorithm would not resolve this, so it is really a "pseudo-wrap". Further, if I add taps to the moving average, it flattens this process out, so I suspect that mathematically, the moving average filter will never reach 0 or 2 pi, though I have never seen an explanation why. Example of an 11-tap:

I find this behaviour fascinating and would be interested in the interpretation of an expert. Does this suggest that features will be distorted at certain "weak spots" in the frequency response? Is it correct to call the phase of a moving average filter "piecewise-linear" rather than linear? I suspect not, given that symmetric FIR filters are analytically shown to have linear phase, but I have a hard time calling this linear.


texts: Is this an infinite response or a finite response?

texts: If you're given $h[n] = 1$ for every even number of $n$, otherwise it'll be 0. Is that considered an infinite response or a finite response and how do you know?


texts: What's the Difference Between LMS and Gradient Descent Adaptation?

texts:
I found algorithms that seems the same to me, but they are described with different names (in field of adaptive filtering).

For example:

LMS - least-mean-squares seems to be GD - stochastic gradient descent
Often the stochastic gradient descent is called just gradient descent what seems to be something different (but still similar) according to wikipedia.
Extension of GD (or LMS) with normalization seems to be always called NLMS, why not NGD?

Why are the used names different? Are the algorithms really the same thing or am I wrong?


texts:
A very simple MATLAB experiment:

f = 200;
fs = 1000;
t = 0: 1/fs : 1;
x = cos(2*pi*f*t);
plot(angle(fftshift(fft(x))));

And here's the output:

Now, made a minor change to the above code snippet; reducing the time duration by just 1 sample, as follows:

f = 200;
fs = 1000;
t = 0: 1/fs : 1 - 1/fs;
x = cos(2*pi*f*t);
plot(angle(fftshift(fft(x))));

And the phase spectrum goes totally crazy:

Questions:

In the first plot, I was hoping to see a zero phase at bin 700 which corresponds to the positive frequency of 200 in this example. That doesn't seem to be the case.
Secondly, I don't understand the linear parts of the graph in plot 1. I do appreciate phase components that could potentially exist due to the so-called numerical noise, but then how could that noise be so 'linear' in phase?
In the second plot, why removing just one sample would have such a drastic impact on the phase plot?
Am I doing something fundamentally wrong here?


texts:
I am looking for the proof that the Hilbert transform of $\displaystyle\frac{\sin(at)}{at}$ is given by

$$\frac{\sin^2(at/2)}{at/2}.$$

How do we prove this? This is a $\operatorname{sinc}(at)$ function. I could not find a proof of this in any textbook so far. Preferably I want to arrive at this result using the analytic signal whose imaginary part would be the Hilbert transform:

$$f_a(t) = 2 \int^\infty_0 FT \Big[ \frac{\sin(at)}{at} \Big] \exp(j2 \pi \nu t) \ d\nu$$

I tried this in a few different ways, but could not get to the answer.


texts:
I was reading a paper about finding optimal morlet wavelet function for CWT, and they were minimizing Shannon entropy.
what is Shannon entropy and it uses?
(specially on continuous wavelet transform)


texts: How to obtain the Spectrogram of a signal

texts:
I am trying to use Matlab to compute the spectrogram of a signal but I can't seem to figure out how Matlab computes the spectrogram.Consider the case where I have a signal x(t) with M samples and I want to compute it's spectrogram using hann windows of length N and using an overlap of K samples.

Because there are M samples and the advance time is (N - K) samples, I would expect the resultant spectrogram to have dimensions of either $N$ by $\frac{W}{(N - K)}$ OR if matlab only provides frequencies up to $\frac{f_s}{2}$ then dimensions to be $\frac{N}{2}$ by $\frac{W}{(N-K)}$ but I am getting weird dimensions.

The code below shows the case where I have computed the spectrogram when M = 32 * 700 = 24400, N = 256 and K = 256 - 32 = 224.

%specgram(data, ChunkSize, SampleRate, Size Of Hanning Window, Overlap (in Samples)
    D = specgram( 1:1:(32*700),256, 8000 ,256 , 224);

produces a 129 * 693 matrix instead of the (256 * 700) or (128 * 700) spectrogram I was expecting.

Result should be an ($M$ by $\text{Number Of STFT Chuncks}$) matrix and
$$
  \text{Number Of Chunks} = \frac{\text{Number Of Samples in Signal}}{\text{Advance Time (in Samples)}}
$$

$$
\text{Number Of Chunks} = \frac{32 * 700}{32} = 700
$$

Have I missed out on something? Any help will be appreciated.


texts: Embedded Image Processing System that sends photos to Smartphone via Bluetooth

texts:
I would like to create an embedded system that takes a Full HD photo by pressing a button,processes and compresses it(for ex. to jpg) and then sends it via Bluetooth to the smartphone - where it's saved in the photo gallery.

I have already started to read "Digital Image processing" by Gonzalez, Woods, and I was thinking if, for my specific purpose, there isn't a better book with which to start.

I am interested to gain understanding about the overall process and the number of components needed for my task and how would they connect with each other.

Do you know any reference material(books/articles/experiments etc.) with which I could start?

I don't believe that I will find everything in one single book, so you can suggest different sources of information for each step needed in the process.

Thanks.


texts: Periodogram and Welch periodogram comparison

texts:
I'm trying to compute the PSD using the periodogram method. I have a signal that is periodical at 5kHz. The frequency resolution sought for is 1.221Hz. It's a PRBS signal that has its energy drop to zero at the generation frequency, i.e. 5kHz in this case. The signals i prety much noise clean.

I tried computing the periodogram and I get what I'm after. Then I tried using Welch's averaging method with N = 500 segment length and overlap of 50% and I get something completely different. I would expect to see an exact copy of the periodogram since there is practicaly no noise. Did I completely misunderstood the method?

What I obtain is this:

I used the following lines of code to estimate the periodogram exactly at the frequencies I want and I scaled the periodogram. ( I have also added .mat files if anyone wants to give it a go)

https://www.dropbox.com/s/pe6qbi883kkocxf/data.mat?dl=0
https://www.dropbox.com/s/2fqmqrg1xpgbibi/freqs.mat?dl=0

load('data.mat');
load('freqs.mat');
1/DT = 5000;
[pxxp,fp] = periodogram(V,rectwin(length(V)),freq,1/DT);
[pxxw,fw] = pwelch(V,rectwin(500),250,freq,1/DT);

semilogx(fp,20*log10(pxxp*8775))
hold on
semilogx(fw,20*log10(pxxw))


texts: Sobel and Prewitt operators decomposition

texts:
I was reading in the wikipedia page of Sobel operator and Prewitt operator that is possible to decompose these two operators (quote form the Formulation paragraph):

  "as the products of an averaging and a differentiation kernel, they compute the gradient with smoothing."

I know that this means that I can rewrite and simplify a Prewitt mask (for example) in this way:

\begin{align*}
\begin{bmatrix}
+1 & 0 & -1 \\
+1 & 0 & -1 \\
+1 & 0 & -1
\end{bmatrix} = \begin{bmatrix}
1\\
1\\
1
\end{bmatrix}
\begin{bmatrix}
+1 & 0 & -1
\end{bmatrix}
\end{align*}

But I don't understand what are in these cases the averaging and the differentiation kernel? And, why it is written that they can compute the gradient with smoothing?


texts: High energy coefficients in 1d DCT for audio

texts:
I've been reading a little bit about how DCT works in image compression using the 2D transform, however I don't know how this applies to audio using a 1D transform.

Suppose I have an audio segment, I mean, some samples of a mono audio vector.

Ai = [-0.0288,-0.0304,-0.0354,-0.0406,-0.0469,-0.0518,-0.0647,-0.0831 ...]

And then I calculate 1D DCT on that vector.

R = dct(Ai)

The result would be something like this:

0.2143   -0.3158   -0.6354   -0.3877    0.5503    0.2330   -0.1690    0.0203

How should I interpret this result? Are the high energy coefficient the ones with a greater value?(I'm still a noob on this). Is the first value the DC coefficient? Is there a special way to read the vector (zig-zag or something)?

Any help will be appreciated, btw I'm working with matlab.
Thanks


texts: PSD of linearly modulated signal using autocorrelation?

texts:
Consider a signal $v(t)$ given by
$$v(t)=\sum_{n=-\infty}^{\infty} b[n]p(t-nT).$$
Assume that $b[n]$ is uncorrelated with zero mean, i.e. $\mathbb{E}[b[n]b^*[m]]=\mathbb{E}[|b[n]|^2]\delta[n-m]$ and $\mathbb{E}[b[n]]=0$.
How can one show,  using the autocorrelation of $v(t)$, that the PSD of $v(t)$ is given by the following?
$$S_{v}(f)=\frac{\mathbb{E}[|b[n]|^2]}{T}|P(f)|^2$$


texts:
I've got a linear time-invariant system $$y[n]=\frac{8}{9}y[n-1]+x[n]$$
which I transformed into a transfer function $$Y(z)=\frac{8}{9}Y(z)*z^{-1}+X(z) =>\frac{Y(z)}{X(z)}=\frac{1}{1-\frac{8}{9}*z^{-1}}=\frac{z}{z-\frac{8}{9}}$$
The pole is positive, which would make the system stabile.
However, I am having some trouble getting the magnitude/phase responses $|H(\omega)|$ and $\theta(\omega)$.
How should I proceed with the system to figure these out? I've read about using the substitution for $z$ as$$Y(z)=\frac{z}{z-\frac{8}{9}}=\frac{e^{iw}}{e^{iw}-\frac{8}{9}}$$
or even continue to transform the exponentials to their $\cos$ + $\sin$ forms, but i am not sure how i should proceed with them.

As an addition, how would the system change if I would to swap $x[n]$ with a sin or a cosine of some value? I assume it would be the same path towards the responses but with different values.


texts: How to extract a smooth contour from a set of points in 3D?

texts:
I am trying to segment a 3D volume. The outcome of all my volume segmentation algorithms is a set of candidate points in 3D space. Now I need to smooth this point cloud and fit a closed surface it. My question is in two parts:

1) How can one smooth a set of rough points to get a set of soft points when demonstrated as a surface?

2) How can I fit a surface to these points to get a closed smooth shape?


texts:
I know I will have a 10% error rate, but this errors are truly random. I was using solomon reed codes, but if a single bit in the code or parity word is bad the whole word becomes useless. If I have eight bytes of data and eight bytes of parity codes then 5 bit error can make my whole message bad. Is there a better algorithm for cases like this?

Edit:

I am trying to send packets of 128bit or 16bytes. 1/2 is data 1/2 is parity. I am having trouble figuring out how LDPC would be more beneficial.

If I have the same 16bytes with 8 bytes of parity and 8 bytes of data how ldpc be better? For example suppose with my example I have 1 bit errors in 5 different bytes. This would make SR unable to decode my data. Would it be the same LDPC or would it be able to recover that. From would I understand from LDPC it would still be unable to decode it.


texts: What is the computational complexity of Spectrum Estimation Algorthms and What they depend upon?

texts:
There are lots of spectrum estimation techniques, each with some pros and cons.
 Algorithms like:

Music
Welch's method
Yule-Walker AR method
Periodogram
modified covariance method
multitaper method (MTM)
Spectrogram using short-time
Fourier transform
Burg method covariance method

and many many other.
 However I do not know which is more proper for embedded system implementation, because I don't know their computational complixy so I am more interested in mathematical estimation of computational complexity of this methods and especially how their complexity is affected with factors like input noise, spectrum occupation, etc?

I'v found Spectrum Estimaton Methods, which greatly explains these methods, but does not compare them from computational complexity point of view.


texts: I am searching for python code or library to detect blank/white images like [png ,jpg ,bmp].


texts:
I have an audio file. I took an extract of the file called variable "CALL", and I took another extract from the file which is just filled with background noise, called "BACKGROUND".
Using matlab, I found the signal to noise ratio of CALL using:
snr(CALL,BACKGROUND); % both variables are of the same length

the signal to noise outputted was     5.3949

I want to add noise to CALL, and lower this signal to noise to 3 for example. How would I go about programming this?

I've only seen examples where people add noise to a perfect signal (i.e. they just create some kind of sinewave, then add noise to that), I have a signal which already has noise, and I want to add more noise, and know its signal to noise value.

Also, is the above method for how I found my first signal to noise value correct? (I'm assuming the background noise is somewhat constant throughout the audio file)


texts: How to estimate for an exponential in WGN

texts:
Question:

Formulating a Maximum likelihood estimator:

So, the likelihood will be

$p(y;\mathbf{h}) = \frac{1}{{(2 \pi \sigma_\eta^2)}^{T/2}} \exp{(-(y - y_0(t))^2)/ 2\sigma_\eta^2}$.

Then, I need to differentiate w.r.t the unknowns and equate to zero.

This is a nonlinear equation in $x$ and cannot be solved directly. Newton-Raphson is a method but it works only for very close initial guess to $x$ and $d$. So, Was thinking how to apply Expectation - Maximization.

I don't know if my approach is correct or not.


texts:
We've studied in signal and systems about odd and even signals. But what is the actual implementation or application of this property?
like for energy and power signals we know If you're using a signal to communicate, then the more energy it has, the farther it will reach (I'm oversimplifying a bit, but in general this is true). The more energy a signal has when it is received, the easier it is for the receiver to recover the information in the signal without (or with few) errors.

  In similar fashion whats with odd and even signal?


texts: How to scale Phase Transition Diagram for Compressed Sensing?

texts:
I want to compute a Phase Transition Diagram as shown here ($A \in \mathbb{R}^{n \times N}$ and $k$ is the sparsity: $\vert \vert x \vert \vert_0 = k $ )

My question is: For $n=1$ I can only compute $k=1$, for $n=2$ I can only compute $k=1,2$ and so on. Therefore the computable number of $\rho$ values is different for each $\delta$ value, but in the plots it is even for $\delta=0$ a continuous, changing plot into the y-direction.
My plots are currently looking like this (yellow means reconstruction with probability 1):


texts: Is there a special case in which wave atoms convert to curvelets?

texts: I'm new to multi scale transformations. I was wondering if there is a special case where wave atom converts to curvelet transform? Can I use wave atom parameters to have curvelet properties?


texts: Normalizing the magnitude of Goertzel Filter c#

texts:
I'm trying to detect some frequencies of 18Khz with sample rate of 44100 [samples/sec].

When I pass a buffer(~4500 samples) to the algorithm it returns a magnitude in a very strange range. Sometimes i get 0.1, or 7 or 250(!) or 154 or 25000 - very strange. I can't decide on a threshold boundary for the filter. Is there any way noramlizing that result? (My result is the magnitude of the complex calculation of real and image)

This is my Goertzel code:

public float goertzel_mag(float[] data,int numSamples, int TARGET_FREQUENCY, int SAMPLING_RATE,int start,int end)
    {
        int k, i;
        float floatnumSamples;
        float omega, sine, cosine, coeff, q0, q1, q2, magnitude, real, imag;

        float scalingFactor = (float)(numSamples / 2.0);

        floatnumSamples = (float)numSamples;
        k = (int)(0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
        omega = (float)(2.0 * Math.PI * k) / floatnumSamples;
        sine = (float)Math.Sin(omega);
        cosine = (float)Math.Cos(omega);
        coeff = (float)2.0 * cosine;
        q0 = 0;
        q1 = 0;
        q2 = 0;

        for (i = start; i < end; i++)
        {
            q0 = coeff * q1 - q2 + data[i];
            q2 = q1;
            q1 = q0;
        }

        // calculate the real and imaginary results
        // scaling appropriately
        real = (q1 - q2 * cosine) / scalingFactor;
        imag = (q2 * sine) / scalingFactor;

        magnitude = (float)Math.Sqrt(real * real + imag * imag);
        return magnitude;
    }


texts: Increase power for particular frequencies

texts:
I am taking an FFT of two signals ( A very clean signal, with little noise, and a very noisy signal).

Signal A

The signal that has very little noise has the following FFT output:

Signal B

The FFT output that has more noise, has the following FFT output:

Is it possible to change the power levels of the frequency bins of Signal B to be more in line of that of Signal A, so that once an iFFT is performed of Signal B, a cleaner version comes out as a result?


texts: Help on understanding Gardner's original paper for QPSK symbol timing recovery?

texts:
I am fairly new to DSP and RF receivers in general, and I was recommended to read Gardner's original paper for symbol timing error recovery for QPSK signals.

I have a question regarding an equation for early-late gates; in the following picture, why is L(r-1) equal to x^2((r-5/4)*T)? It seems to me that it should be x^2((r-3/4)*T). Am I missing something?

Also, intuitively, if the spacing for the Early and Late is 1/4 each, then the spacing between Early and Late should be 2 * 1/4 T = 1/2 T, not 1 * T.


texts: Fourier Transform of triangle function $x(t)=\Delta\left(\frac{t-1}{2}\right)$

texts:
Can you please tell me if my working is right for the Fourier Transform of this function: $$x(t)=\Delta\left(\frac{t-1}{2}\right)$$

My workings are:

I have used the fourier transform standard results. Please could you tell me if I which I have to apply:

First the $F(t-1)$, then $F(t/2)$
First the $F(t/2)$, then $F(t-1)$
First the $F(t/2)$, then $F(-1/2)$

Please tell me where I have gone wrong.

My final answer is:

$$x(t)=\mathrm{sinc}^2(\omega/2)e^{-j\pi f}$$


texts: Use of pseudo-random binary sequence (Gold, Kasami, etc.) in frequency hopped code-division multiple access (FH/CDMA)

texts:
I have a clear idea of how a PN-sequence is used in a Direct Sequence Code-Division Multiple Access(DS-CDMA) system with BPSK. In such a system, a different pn-sequence is assigned to each user and each transmitted bit is XORed with the binary pn-sequence.

Now my question is:

For each symbol time, how are this binary pn-sequences used in a FH/CDMA system to chose which hop frequency to use if there are more than two possible hop frequencies?

If I have 2 possible hop frequencies to use, the answer is trivial, because each bit of the pseudo-random sequence would tell which of the two hop frequencies must be used. But when having more than two frequencies... How is this selection done?


texts:
One of my tutorial questions for communication systems asks me to find the time function $x(t)$ which has the Continuous-Time Fourier Transform: $$X(\omega) = \frac{3}{(1+j\omega)(2-j\omega)}$$ So far all I have been able to do is rationalize the function so by multiplying the numerator and denominator by $(1-j\omega)(2+j\omega)$ and then substituting the resulting function $X(j\omega)$ into the formula for the inverse Fourier Transform. I can't find examples like this anywhere and so I can't evaluate where the upper and lower limits should be in order to solve this.

Does anyone have some knowledge to shed some light on this for me.

Simon.


texts: discrete fourier transform in Matlab - theoretical confusion

texts:
I have a periodic term $V(x) = \sum_K  \exp(iKx) V_K$ where $K =2\pi n/a$ where $a$ is the periodicity of the term and $n =0,1,2,3....$

Now I want to find the Fourier coefficient $V_K$ corresponding to a particular $K$. Suppose I have a vector for $V(x)$ having $10000$ points for $x = 0,0.01a,0.02a,...a,1.01a,....2a....100a$ such that the size of my lattice is $100a$.

FFT on this vector gives $10000$ Fourier coefficients. The $K$ values corresponding to these Fourier coefficients are $2\pi n/(10000*0.01a)$ with $n=0,1,2,3,...9999$. But my $K$ had the form $2\pi n/a$ due to the periodicity of the lattice.

What am I missing ?


texts: Kalman Filter: Why do we decrease the state uncertainty regardless of the current measurement?

texts:
I'm struggling with fully understanding the concept behind a Kalman filter.
For the sake of simplicity, let's ingore the input variable $u$ and assume constant process $Q$ and measurement noise $R$.
We can now predict the new state by using the state transition matrix $F$ and the old state estimation
$$x^-_k = Fx_{k-1}$$. Same does apply for the estimate uncertainty matrix : $P^-_k = FP_{k-1}F +Q$. Now we observe some data and depending on the ratio between the estimate uncertainty and the measurement uncertainty, we update our belief of the state estimation by calculating the innovation
$$x_k = x^-_k + K(z_k - Hx^-_k)$$
where $z_k$ is the current measurement, $H$ is the  observation matrix and $K$ is the kalman gain $P^-_kH^T(HP^-_kH^T+R)^{-1}$.
To incorporate the new measurement the state uncertainty is updated with $P_k = (1-KH)P^-_k$.

Do I see it correctly that $P_k$ decreases with each iteration because $ (1−KH)≤1$? Assuming that this is true, doesn't this mean that the uncertainty in the state estimation becomes very small with enough time steps and eventually tends towards $Q$ since the values within KH become smaller with each step?

If we are at a point where the values in $P$ are very small, which yes means that the uncertainty in our estimated state and the estimated measurement is very small, then the influence of the measurements is probably very small as the Kalman gain becomes small - I.e. even if the innovation $z_k-Hx^-_k$ is very large, we trust the prediction of the model, although perhaps the true observation would have been correct?

Suppose the values in $P_k$ are very small and we get measurements that are very different from the previous ones, but (for whatever reason) are true. The predictions of the kalman filter will deviate strongly from this and still have a small variance, simply because according to the previous data the prediction of the KL is more probable than the measurement? I.e. At a certain point  we have only predictions with a small variance, no matter whether they correspond to the truth and I see now way how the values in $P_k$ could grow again. But I think they should grow because if the innovation is big for some iterations, then there must be something wrong with the predictions

Thankful for any tip or reference


texts: Possible to break out complex WiFi signal?

texts:
Is it possible to break out the i and q of a complex WiFi signal?  If I have a laptop with a standard WiFi card in it, could I break out the i and q in analog, or upon sampling in an FPGA via digital means?

It must be possible since it is being done on a WiFi receiver, right?


texts: Constructing frequency filter from a 2d spatial filter

texts:
Given a 2d spatial kernel $h$, I want to construct the corresponding frequency filter $H$. In theory, $H$ is simply obtained by calculating DFT of $h$ using a FFT algorithm. But in practice, there are some considerations to be considered. First, if $H$ is used for filtering, zero-padding is required. Second, after zero padding, the original center of $h$ (i.e. center of $h$ before zero padding) must circularly-shift to the upper-left corner.
I know that zero padding is used to overcome wraparound errors. And I have observed, by experiment, shifting the original center to the upper-left corner prevents shifting the filtered image. But my question is why this is the case? That is, why not shifting the original center, leads to shifting the final image; and how shifting the original center resolves the problem?
Thanks


texts:
I have implemented Haar Function in MATLAB based on formulas stated here. However this formulas generate functions only between 0 and 1 (on [0..1] interval). I need to generate such functions periodically from -infinity to +infinity for instance with period = 1. It should depend on length of my time vector. For example if time vector on interval [0..10] I should generate Haar function vector on the same interval [0..10]. I hope you get a point. Here is my MATLAB code.

function [ f ] = haar(x)
%Haar function
%   x: time vector
vectorLength = length(x);
f = zeros(1, vectorLength);
for k = 1:vectorLength
    if (x(k) >= 0 && x(k) < 1/2)
        f(k) = 1;
    elseif (x(k) > 1/2 && x(k) <= 1)
        f(k) = -1;
    end;
end;
end


texts: SSD object detection - Class score and bbox coordinates computing

texts:
How does SSD ( single shot multibox )object detection model computes the the confidence score and the bounding box coordinates, referring to SSD paper, it said that this process is done with the aid of small convolutional Kernels. The kernels will output four real values (bbox offsets) and a vector with class scores. Can someone explain how this process is done with more clarification?

SSD: Single Shot MultiBox Detector available at:

https://arxiv.org/pdf/1512.02325.pdf

Thanks in advance.


texts: Multidimensional FFT with mixed directions

texts:
Let's say I have three-dimensional complex array $X[m,n,p]$, $m \in \{0, 1, \ldots, M-1\}$, $n \in \{0, 1, \ldots, N-1\}$, $p \in \{0, 1, \ldots, P-1\}$; i.e. $X \in \mathbb C^{M \times N \times P}$.   I would like to take a three-dimensional DFT of this array; along dimensions $m$ and $n$ I would like to use the forward transform, and along dimension $p$ I would like to use the inverse transform.

Every multidimensional FFT implementation I've looked at (FFTW, MKL, numpy.fft) just has one parameter to specify as "forward" or "backwards/inverse", implying that the same transform is applied in each dimension.

Is there an "easy" way to do this mixed-direction FFT directly without having to break it up into two steps (i.e. doing a 2D forward transform and then a set of 1D inverse transforms)?

Either a conceptual answer (e.g. "conjugate this then do this") or a "programming" answer (e.g. "use this flag in MKL") would be fine.


texts: How to calculate DTFT of cosine function divided by n

texts:
I'm having a hard time to calculate the next function, and I don't really know Matlab good enough to calculate it there.
Help would be appreciated:
$$h[n]=\frac{A_1 \cos[\theta_1(n-N/2)]}{n-N/2}$$


texts: Radar signals features for machine learning

texts: I have project in which I have to classify materials based on the spectra of received  radar signal. I want to create feature vectors of the specific materials. Can anyone share his ideas. ?


texts: Expressing 2N point DFT in terms of N point DFT

texts:
I have a problem with expressing odd samples of X2 in terms of X1. I understand that the resulting DFT will be more precise in terms of expressing the exact spectrum of signal x[n], due to more samples. Moreover I know that even samples of X2(k) are copy of the spectrum of X1(k), however I do not know how to mathematically compute the odd ones.

The task is as follows:

Let $x[n]$  be  a  periodic  sequence  with  period $N_1$ .   Thus $x[n]$  is  also  periodic  for  period $N_2 =  2 N_1$ .   We  may compute $X_1[k]$:  $N$-point DFT of $x[n]$ and $X_2[k]$:  $2N$-point DFT of $x[n]$.

Express $X_2$in terms of $X_1$ Hint: it is easy with even samples $X_2[2m]$ , harder with odd ones $X2_[2m+ 1]$ , for $m$ integer.


texts: How to perform median filtering in multi-label prediction

texts:
I'm trying to perform automatic pitch tracking for piano. It's easy to detect the note when there's only a single note pressed. There are 88 keys on the piano and I set the result to one of them. There may be misclassifications due to octave shift. For example, instead of 40, I find 52. They are both C, but 52 is one octave up. I can remove most of these octave shifts by performing median filtering on the results.

But it is confusing when I move to polyphonic pitch tracking. When there are more than a single note pressed together, there can be varying number of pitches detected at time t.

For example, let labels(f,t) be 1 if note f is detected at time t and let the binary matrix labels it be as follows

labels = 0 0 0 0 1 ...
         1 1 0 1 0 ...
         0 1 0 0 1 ...
         .............
         0 0 1 0 1 ...
         .............
         0 0 0 0 1 ...

Let the detected notes be as follows correspondingly:

labelsAsNotes = 2; 2 3; 14; 2; 1 3 14 88; ...

How can I do median filtering in this case? Do you have any other recommendations to remove octave shifts?


texts: What is the relationship between group delay and propagation delay for a filter?

texts:
I was trying to design a filter for a system and I was able to obtain the group delay plot in matlab using grpdelay. The group delay output is in samples. However, I was not able to relate that data to the propagation delay of the signal through the filter. Is there a formula between these two quantities?
Also, if I were to change the magnitude of the input signal in a step (keeping the frequency constant), how long would this change take to propagate to the output? Would that depend on the group delay or phase delay or would it depend on the filter order?
For example, if the filter cutoff is 5000Hz and I have unity gain at 500Hz and Group delay of 15 samples and Phase delay of pi radians, then if I send a 500Hz signal through the filter, after how many samples would I get it at the output? Assuming sampling frequency is 10000Hz.


texts: Spectral density of signal written as quadratures

texts:
Suppose I have a real-valued random signal $f(t)$ written as
$$f(t) = I(t)\cos(\Omega t) - Q(t)\sin(\Omega t)\, .$$
What is the relationship between the spectral densities of $I$ and $Q$ and the spectral density of $f$?

One approach is to use the fact that
$$\int_0^\infty \frac{d\Omega}{2\pi} S_f(\Omega) \cos(\Omega \tau)=
\langle f(t)f(t+\tau)\rangle_t \, .$$
But so far I haven't understood how to make progress.


texts: How sampling aperiodic signal will result in periodic repetitions of the same

texts:
I am reading "Digital Signal Processing" - Proakis and often read that sampling guarantees periodicity (not exact as read)

But I wonder how sampling aperiodic signal will result in periodic repetitions of the aperiodic signal?


texts: Understanding the Windowing Method in PSD Calculation

texts:
I have some issues understanding the use of a window (it doesn't matter which one) to calculate the Power Spectral Density of a signal.

For example say I want to use the hann window, with $N = 1024$, knowing that my signal length is $X > N$, $X$ is arbitrary.

How should I calculate the PSD :

Take 1024 samples from the signal, multiply them in the time domain  with window function. Then transform those 1024 time values. After that take the next 1024 samples of the signal and do the same thing, then add the result to the spectrum that I got the first time and so on?

That's the way I understood this, but to be honest I'm not sure that it is a right way?


texts:
What are Leslie speaker / "rotary" style of digital audio processors doing?

I think a simple one would model:

Delay between speaker and a pair (stereo) of microphones.
The tonal change when the speaker rotates (LFO modulated filter and amplitude modulation).

But what (more) is there really?


texts: vibration analysis in spherical coordinates

texts:
I'm trying to understand some data from a vibration sensor, in units of milli-g and in cartesian coordinates.  One idea I had was to analyze the data in something like spherical coordinates, to understand if most of the vibration is radial, or how much is transverse along the angle dimensions.  Can the FFT algorithm be used on data converted to spherical coordinates?  Or are there similar algorithms for converting time series spherical coordinates to frequency domain?  Can the FFT be used on the radial coordinate?

I have a physics background, but its been awhile.  I'm trying to remember somethings from old courses and wanted to check if this is valid.


texts:
I'm trying to implement a QPSK burst modulation and demodulation through a channel and I'm now considering the synchronization part.

I've read that Barker sequence were excellent option for this as their autocorellation looks like a dirac.

So far I know that they're composed of 1 and -1 and that the longest found is 13 characters long.

My question is the following:

Is there some adaptation of a Barker sequence for QPSK to make use of
the complex representation or do you simply send a classic sequence
using only 1 and -1 while you could use also i and -i?


texts: A little bit better than zeroing bins in order to filter

texts:
The topic of zeroing frequency bins before IFFT has been discussed here : Why is it a bad idea to filter by zeroing out FFT bins?
It was very helpful to understand why it is a bad idea.

What do you think about, instead of zeroing, multiplying by a window (in the frequency domain!) which is 0 at the centre and smoothly goes to 1 outside of a small interval ?

t = linspace(0, 1, 256, endpoint=False)
x = sin(2 * pi * 3 * t) + cos(2 * pi * 100 * t)
X = fft(x)
window = 1 - scipy.hamming(128)
X[64:192] *= window
y = ifft(X)

Then it would not be a "rectangular" zeroing, but a "smoothed zeroing" in the frequency domain. What do you think about such a filter ?


texts: Matlab : Exploring the effects of co-efficient quantization

texts:
b = [0.9355,-3.6135,5.3603,-3.6135,0.9355]
a = [1,-3.7339,5.3562,-3.4929,0.8752]

Use MATLAB to quantize the coefficients in the vectors a and b to 16 bits of precision. First, take each vector and divide by the smallest power of two such that, after division, the magnitude of the largest filter coefficient is no larger than 1. Next, multiply the resulting vectors by 215, and round to the nearest integer. (You can use the round function.) Divide the resulting vectors by 215. The resulting values should be no less than -1 and no larger than 1. Finally, multiply these values by the power of two by which you originally divided.

I do not understand what that smallest power of two is! Can someone help me out with the code?
Thanks!


texts: How to construct a preemptive digital filter to neutralize an unwanted filter

texts:
I'm trying to think of how to construct a filter based on the following scenario. I have some time trace with a certain known power spectral density (I can verify this with a periodogram of the time trace). The time trace then goes into a machine that realizes it, which does so with some level of imperfection. Treating the device as a black box, it is essentially a filter that is applied to the time trace, and produces a slightly adjusted power spectral density. This power spectral density I can then measure with a spectrum analyzer.

So I'm wondering, can I apply a preemptive filter to my time trace, before it goes into the black box, so that I get my desired power spectral density? I suppose I should have all the information I need as I know the power spectral density before and after the black box, which should somehow allow me to construct a type of inverse filter. How exactly I am not sure of however, and that is my main question.

I should also note that I tried a much less clean method. I took the power spectral density after the black box, and simply mirrored it, and used this as a preemptive filter. This actually seems to work quite well in the sense that the resulting PSD has the correct shape (it was sloping downwards before while it should be straight, which this method achieves), but the amplitudes are no longer correct and rescaling the data is problematic. So I figured a more analytic approach would be required.


texts: Ideas on Matrix Factorization / Transformations for $ {L}_{1} $ Minimization

texts:
I am starting with a typical $\ell_1$ basis pursuit problem:

$$
\min_{\mathbf{x}} \Vert \mathbf{x} \Vert_1 \quad \mathrm{s.t.} \quad \Vert \mathbf{ERx} - \mathbf{y} \Vert_2 \leq \epsilon,
$$

where $\mathbf{R}\in\mathbb{C}^{M \times P}$, and $\mathbf{E}$ is an $M' \times M$ row selection matrix that randomly chooses $M'<P$ rows from $\mathbf{R}$. It is also known that $\operatorname{rank}(\mathbf{ER})=M'$, so that the constraint system is full-rank but under-determined.

Since both $M$ and $P$ are on the order of $10^7$, the above problem is impractical because of memory and computational limitations. An algorithm exists that will let me efficiently compute $\mathbf{R^{\mathrm{H}}}\mathbf{E^{\mathrm{H}}}\mathbf{y}$, so I would like to see if I can combine this with a sparse approximation (described below) to solve the following equivalent problem:

$$
\min_{\mathbf{x}} \Vert \mathbf{x} \Vert_1 \quad \mathrm{s.t.} \quad \Vert \mathbf{R^{\mathrm{H}}}\mathbf{E^{\mathrm{H}}}\mathbf{ERx} - \mathbf{R^{\mathrm{H}}}\mathbf{E^{\mathrm{H}}}\mathbf{y} \Vert_2 \leq \epsilon,
$$

Since $\mathbf{E}$ is a row-selection matrix, the product $\mathbf{E^{\mathrm{H}}}\mathbf{E}$ is a diagonal matrix with a $1$ at element $(i,i)$ if row $i$ was selected and $0$ otherwise. Furthermore, given the nature of $\mathbf{R}$, even though the product $\mathbf{R^{\mathrm{H}}}\mathbf{E^{\mathrm{H}}}\mathbf{ER}$ is rather large at $P \times P$, it is also approximately sparse with only a relatively small number of terms with significant modulus per row/column. This makes me think that the above form is doable since I could compute just the sparse approximation and avoid the unrealistic memory and computation requirements.

However, $\mathbf{R^{\mathrm{H}}}\mathbf{E^{\mathrm{H}}}\mathbf{ER}$ has size $P \times P$ but is only of rank $M'$ and this is causing problems with the  basis pursuit algorithms as they expect the rows of the constraint matrix to be
linearly independent. Hence, my solicitation for input.

Are there any factorizations and/or transformations that could turn the constraints of the second form of the problem into under-determined yet full-rank system? Please do keep in mind the large sizes of these matrices. Any other suggestions/comments are also welcome. We can assume that $\mathbf{R}$ has rank $M$ or $P$, as there is freedom in its design.


texts: Magnitude of Power Spectral Density with Different Sampling Frequency

texts:
For exponential signals (sine or cosine), if the sampling frequency $f_s$ is equal to the length of signal, $N$, the magnitude of psd for each sine signal is proportional to amplitude in time domain, while if $f_s$ is not equal to length of signal $N$, for example here, $f_s=1000 \text{Hz}$, the magnitude of sine signals in psd (frequency domain) is not proportional to their original signal. The MATLAB code is provided.

Can you tell me the reason? Is this something related with leakage?

N=256;
x0=zeros(N,1);
fs=1000;
idx = (0:(N-1))'; % Indices
f=20*[1 2 3 4];
A=[1 2 3 4];
K=4;
for ii=1:K
    x0=x0+A(ii)*exp(sqrt(-1)*2*pi*f(ii)*idx/fs);
end
X=fft(x0);

figure;
plot(abs(X).^2)


texts: Should you include division by coefficient of variation in WV reconstruction for its normed 1D signal?

texts:
I make a new ECG time series from WV spectrum of original signal and its L2 energy normalisation.
I am thinking if the reconstruction step benefits from covariance at each time point.
I take later a norm of the matrix, which creates the new time series.

result there is no benefit of adding cov postop into the L2 new signal - verified, strong evidence
now, testing if the application of cov on the initial energy spectrum helps

Nonlinearity

I think the coefficient of variation should not be applied to the initial spectrum because there exists a bijection between the two data sets, although a nonlinear one.

  The problem here is that neither sensor is better in this case, because the data sets are direct conversions of each other, but the CV of each data set is different: 0.79 versus 0.42.

Hovever, I am not sure if the nonlineary of the problem case makes the thing different.

Source: https://en.wikipedia.org/wiki/Coefficient_of_variation


texts: IIR Filter structure implementation dynamic range

texts:
I want to implement an IIR filter.  The filter is going to implement by cascading second order sections. Below there are 3 pictures that implement this filter using:

Denominator - Numerator using Tapped Delay Line Structure 2:
Numerator -Denominator using Time Delay and Accumulate Structure 3:
Denominator - Nominator using Time Delay and Accumulate Structure

(1)

(2)

(3)

a,b: coefficients

I want to examine the efficiency of the structure due to dynamic range and ease of implementation. I also want to use only one multiplier in the architecture. I think it is easier to implement the structure with less delayers but what is the dynamic range and how i can find it??


texts: Techniques Used To Process Accelerometer Recordings

texts:
I'm currently working on a project which involves using Accelerometer data. Having searched the internet for techniques in processing the data, I am unable to find simple techniques. Is there any suggestions, of what techniques to use in processing recordings? The data recorded is to do with body motion, walking up-stairs, down-stairs, falling etc.

The advantages and disadvantages of techniques to processing the data would also be welcome.


texts: Difference between Fourier-Transform and FFT of rectangular pulse

texts:
I'm trying to find a link between the Fourier-Transformation of aperiodic Signals and the FFT of them.
So to start with a basic example, let's take a rectangular pulse with width 0.1s and amplitude of 1 shifted by 0.05.
Using correspondency, i can calculate the expected spectrum:
$X(f) = 0.1 \cdot sinc(0.1f) \cdot e^{j 2 \pi f \cdot 0.05} $

But now, when i generate the Signal with the follwing Matlab-Code:

f_abt = 50e3;
x=0:1/f_abt:1;
y=zeros(1,length(x));
for ii=1:length(x)
    if x(ii)<=.1
        y(ii)=1;
    end
end

And calculate the spectrum of it, the result is depending on the length of the signal. So, when i compute the one-sided spectrum from the signal generated above (1s duration), i get:

Then, when i put the signal length to 2s (everything else unchanged):

x=0:1/f_abt:2;

I get the following spectrum:

I guess the difference comes from the FFT-Algorithm i use. When doing FFT, i normalize the Values by Nfft, so it makes complete sense that my amplitudes change when i change the Signal length.

My question is: How do i get the right spectrum and how do i know it is right, e.g. when i can't calculate it 'by hand' using correspondencies? I'm having issues finding the connection between my "real", time limited signal and its FFT and the "theoretical" rectangular pulse.

Code I use for calculation of one-sided spectrum:

function [f_xa, mag, phase] = calc_fft_f(ta, xa)
N_a = numel(xa);
fft_xa = fft(xa);
P2_norm = fft_xa/(N_a);
if (mod(N_a,2))
    P1_norm_single = P2_norm(1:ceil(end/2));
    P1_norm_single(2:end) = 2*P1_norm_single(2:end);
else
    P1_norm_single = P2_norm(1:(end/2)+1);
    P1_norm_single(2:end-1) = 2*P1_norm_single(2:end-1);
end

mag = abs(P1_norm_single);
phase = rad2deg(angle(P1_norm_single));
Fsa = 1/(ta(2)-ta(1));
f_xa = Fsa*(0:(length(mag)-1))/N_a;

end

Thanks in advance!


texts:
Trying to modulate audio with OFDM data.

My fs = 44.1khz, i'm using 64 sub-carriers, but i'm only using 5 sub-carriers for data, so my bandwidth is about 3.445khz. I'm not sure how to properly write the all real data to a wav file and then reconstruct the complex signal when I read it back. Can someone point me where I'm messing up?

Matlab code:

 close all
clear all

%fft length
N=64;
NN=1024;

%carrier frequency
fc = 15000;

%create the data to send
qpsk1 = (floor(2*rand(1,2))-.5)/.5 + 1j*(floor(2*rand(1,2))-.5)/.5;
qpsk2 = (floor(2*rand(1,3))-.5)/.5 + 1j*(floor(2*rand(1,3))-.5)/.5;

inputiFFT = [zeros(1,29) qpsk1, 0, qpsk2 zeros(1,29)];
outputiFFT = ifft(inputiFFT,N);

%spin data to carrier
y3 = outputiFFT.*exp(1j*2*pi*15000/44100*(0:length(outputiFFT)-1));
%despin to baseband for sanity check
y4 = y3.*exp(-1j*2*pi*15000/44100*(0:length(y3)-1));

figure
subplot(211)
hold on
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(outputiFFT,NN)))))
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y3,NN)))))
hold off
title('Data - blue, Data at carrier - r')
subplot(212)
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y4,NN)))))
title('sanity check down mix')

%upconvert to 15khz and write to wav (all 'real')
towav = real(outputiFFT).*cos(-2*pi*fc/44100*(0:length(outputiFFT)-1)) + imag(outputiFFT).*sin(-2*pi*fc/44100*(0:length(outputiFFT)-1));

%write to wav then read from wave
wavwrite(towav, 44100, 'ofdmTest.wav');
[y_wav,FS] = wavread('ofdmTest.wav');

figure
hold on
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(towav,NN)))),'color','b')
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y3,NN)))),'color','r')
hold off
title('Theoretical data up converted to carrier -r , data written to wav - b')

%convert to I-Q
yRx = (y_wav.*cos(-2*pi*fc/44100*(0:length(outputiFFT)-1)).'+ 1j*y_wav.*sin(-2*pi*fc/44100*(0:length(outputiFFT)-1)).');

figure
hold on
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(yRx.',NN)))),'color','b')
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y_wav,NN)))),'color','r')
plot((-0.5:1/NN:.5-1/NN)*44100, ((abs(fft(y3,NN)))),'color','g')

hold off
title('yRx -b, data written to wav -r, theoretical - g')


texts: How left out down-sampling in 3D-DWT just on the Z direction

texts:
I am working on image processing by a  3D wavelet transform.

I have a problem with the classification size of the wavelet coefficients. As you know, when we apply 3D_DWT on the image (for example: Hyperion of the area by 256*256*128 size) the sub-band size would be decreased (for LLL in 1l level the size is 128*128*64).

How I could apply down sampling just for the Z direction (spectral direction)?

It means this size of the sub-band should be 256*256*62 after one level wavelet applying. I want to do it in the Matlab script (dwt3.m & wavedec3.m), but when I removed the down sampling, it was applied in all three directions.


texts: Meaning of Q and I in all-software QPSK communication systems?

texts:
Q and I seem to come up a lot in discussions of digital communications, such as the comprehensive answer provided to this question:

FFT window phase syncronising

It seems like the Q and I channels ultimately are combined into a sin wave of a given frequency, phase and amplitude.  If the communications system is all software (for example in a software-defined radio) is it necessary to worry about "Q" and "I"?  Is this just an artifact of when receivers and transmitters were electronics based?

To make the question more concrete, the link provided talks about putting the even bits onto Q and the odd bits onto I.  If this is the same thing as creating a sin wave of phase 0, 90, 180, or 270, then why bother with this seemingly unnecessary notation?  Perhaps the pulse shaping (like an RRC filter) needs to be independently applied to the Q and I channels?


texts: Deconvolve a FIR-filter using limits? Or: Are limits distributive with the inverse fourier transform?

texts:
Given $y[n] = f[n] * h[n]$ where $y[n]$ and $f[n]$ are two one-dimensional discrete  signals that are given, find out:

if $h[n]$ is a FIR-filter
if it is a FIR-filter, what its kernel is.

Obviously, this can be done by doing deconvolution: $h[n] = y[n] *^{-1} f[n]$.
Performing deconvolution (written here as $*^{-1}$) in the time domain is of course computationally expensive. However, it is possible to go to the frequency-domain instead:

$$h[n] = y[n] / f[n] = \mathcal{F^{-1}}(Y ./ F) $$
(Where $ ./$ is pairwise division and $Y$ and $F$ are the (discrete) fourier transformations of $y[n]$ and $f[n]$ respectively, that is $Y = \mathcal{F}(y[n]$), $F = \mathcal{F}(f[n])$).

Many existing algorithms use this property, since pairwise division has a lower time-complexity) than direct deconvolution.

However, these algorithms claim to struggle with zero-valued frequencies: If $F$ contains a zero, then we will divide by zero during the pairwise division.

My question is: Is it possible to calculate the outcomes using limits here? After all:

\begin{align}
\lim_{k\to 0^+} \frac{x}{k} &= \lim_{l\to \infty} l \tag{$0^+$}\\
\lim_{k\to 0^-} \frac{x}{k} &= \lim_{l\to -\infty} l  \tag{$0^-$}\\
\mathcal{F}(x) &= \int_{-\infty}^{\infty} f(x)\ e^{-2\pi i x}\,dx\tag{FT}\\
\mathcal{F^{-1}}(x) &= \int_{-\infty}^{\infty} f(x)\ e^{2\pi i x}\,dx\tag{IFT}\\
\int_{-\infty}^{\infty} f(x)\,dx &= \lim_{h \to 0} \sum{\frac{f(x)+f(x-h)}{2}} \tag{sum limit integral}\\
\end{align}

So the fourier transform (and the inverse fourier transform) are defined in terms of integrals, which are themselves defined in terms of limits.

Limits distribute over summation, so it seems to me that it would be possible to simply use IEEE 754 floating-point numbers (which follow the limit-rules for division by zero) to find a correct answer, which would allow for a very fast computer-implementation of a deconvolution algorithm.

As an aside, I wonder if the fact that $h[n]$ is a FIR-filter (and linear, time-invariant and causal) is useful in some way.

Is my reasoning correct, or is it flawed somewhere?


texts: After a BPSK modulation we have the upsampling and the pulse shaping filter. But for a FSK, do we have upsampling and pulse shaping filter, or just the pulse shaping filter?


texts:
So I have a digital linear chirp signal:

$$ x[n] = \sqrt{\frac{2}{N}} \cdot \cos\Big( 2 \pi \frac{f_{0}}{f_{s}}n + \pi \frac{B}{N} n^{2}\Big) $$

Where $f_{0}$ is the initial frequency, $f_{s}$ is the sampling frequency, $B$ is the bandwidth and $N$ is the number of samples.

I want to simulate a Doppler shift. I am not sure exactly how to do this.

I know I need to multiply by initial frequency $f_{0}$ by the Doppler shift $d$, and do the same with the bandwidth. This will give me:

$$ x_{d}[n] = \sqrt{\frac{2}{N}} \cdot \cos\Big( 2 \pi \frac{d \cdot f_{0}}{f_{s}}n + \pi \frac{d \cdot B}{N} n^{2}\Big) $$

I know Doppler compresses or expands my waveform in time, and I know how to simulate the Doppler shift in the case of a simple sinusoid. However, I am a bit confused here.

Do I need to change the number of samples $N$ as well? If so, by how much exactly?

Thanks


texts: Reconstruction of audio signal from Spectrogram: Implementation

texts:
Is there somewhere an implementation available for the reconstruction of audio from a spectrogram? (e.g. based on this approach).

I could not find anything, but I think there should be  an implementation somewhere. Where can I find it?

Thank you


texts: Seemingly no division operation is done of the diagram of the transfer function

texts:
Given the below system.
$$  y\left[n   \right] =2 x\left[n \right] -x\left[n -1\right] +0.5 y\left[n -2  \right]  $$
We'll find out the transfer function of it.
So, easiely using z-transform,
$$  Y(z) =2 X(z) -X(z) z^{-1} +0.5 Y(z) z^{-2}  $$
$$ Y(z) \left\{ 1-0.5 z^{-2}  \right\} =X(z) \left\{ 2-z^{-1}  \right\}   $$
$$ H\left( z  \right) = \frac{Y(z) }{X(z) }= \frac{2-z^{-1} }{1-0.5 z^{-2} }   $$
gained it.
Currently I can't get the below diagram.
I can't see component which operates division so it is weird that the diagram is applied to the transfer function.
About the blue domain of below, currently I can't get how the 1 of the denominator of the transfer function is attained.

What I've been missing or what should I study for next?


texts: Phase of matched filter for chirp signals

texts:
I am using the following system: Ultrasonic chirp signals with central frequency 100kHz are being transmitted using ultrasonic sensor. After the transmission through the medium(air/solid/liquid/combination) they are being received with hydrophone. The system is being controled by DSP and signal generation and acquisition is obtained via DAC and ADC. After reception we are performing demodulation and then complex matched filtering. Our aim is to detect the time of arrival between transmission and reception. So far I have been only observing the peak of absolute value of matched filter output as the moment that corresponds to time of arrival, but if I am correct the phase of the matched filter output can also be used.
My question is, should the phase of matched filter output be zero at the moment of the peak in the absolute value of matched filter output?


texts: Detecting lines with houghlines in matlab

texts:
I have This image which is basically made of sloped lines (almost horizontal and parallel)

I tried exactly the code in matlab documentation to detect those lines (except the rotation part) but I got a very strange result (vertical lines detected and no single horizontal line!):

I really can't find explanation of this output. I tried the same code on a different image and it was ok. I tried to change all the possible parameters in houghlines and houghpeaks but I always get vertical lines even though my image has no vertical lines or anything similar to be vertical!

I need to detect the parallel lines in the image or somehow get their slope. If you think there's something wrong in my image or additional modification I should do to the example code. or if you have any better idea that would be great!


texts: Example of an LTI system with complex impulse response

texts:
In general, I know that the impulse response $h(t)$ of an LTI system can be complex. However, all of physically realizable, useful systems I've come across have purely real impulse responses. I did a web search in an attempt to find a useful, non-pathological system with a complex impulse response, but was not immediately successful.

Could anyone provide an example of a system with a complex impulse response, and state what applications it has? Links to such examples would also be appreciated.


texts:
What kind of FSK signal is this and what demodulation technique can I use to demodulate it. Why is the phase of the 1's changing like that and how does that effect how it needs to be demodulated. I will be programming the demodulation in python.

I am trying to low pass at 900Hz and then generate the envelope but I am not getting the desired result.

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal
from scipy.fftpack import fft, rfft, rfftfreq, irfft
import scipy.signal.signaltools as sigtool

from scipy.io import wavfile
fs, data = wavfile.read('jx3p.wav')

h = signal.firwin(numtaps = 300, cutoff = 900, fs=fs)
data = signal.lfilter(h, 1, data)

data = np.abs(sigtool.hilbert(data))

import matplotlib.pyplot as plt
plt.plot(data)
plt.show()

UPDATE

Here is the code to multiple the signal by the same signal with a 12 sample delay and then add the resulting signal with a 3 sample delay of itself.

import matplotlib.pyplot as plt
from scipy.fftpack import rfft, rfftfreq
from scipy.io import wavfile
from scipy.signal import blackman

fs, data = wavfile.read('jx3p.wav')

x = 286600
y = x+22050
data = data[x:y]

m = max(abs(data))
data = [d/m for d in data]

delay = 12
data1 = tuple(a*b for a,b in zip(data[delay:],data))

filter_delay = int(delay/4)
data2 = tuple(a+b for a,b in zip(data1[filter_delay:],data1))

plt.plot(data[:300],"r")
plt.plot(data1[:300],"b")
plt.plot(data2[:300],"g")

plt.show()


texts: Generating a time series given a transfer function

texts:
I'm trying to work my way through a paper I found online, titled "Three Models of Wind-Gust Disturbances for the Analysis of Antenna Pointing Accuracy" by W. Gawronski, 2002.

http://ipnpr.jpl.nasa.gov/progress_report/42-149/149A.pdf

In the Appendix of the paper, Section I, the author outlines an approach to generate a time series.

  (1) Apply white noise to the ﬁlter, H, with the following transfer
  function: $$H_d = \dfrac{0.1584 \cdot z^3 −0.3765 \cdot z^2 + 0.2716 \cdot z −0.0534}{ z^4 −2.9951 \cdot z^3 +3.0893 \cdot z^2 −1.1930 \cdot z + 0.0988}$$ and obtain time series of wind $\Delta v(t)$.

How do you apply white noise to a filter? How do you use the transfer function to create a time series?


texts: Understanding CQT (Constant-Q Transformation) parameters for piano AMR (automatic music recognition) with Librosa Library

texts:
I am doing an automatic music recognition project with a deep learning model. For my data preprocessing, I am trying to calculate the Constant Q Transform for Polyphonic 88-key piano audio using Python's Librosa library. However, I do not understand what I should set fmin, n_bins, bins_per_octave to in Librosa's CQT method to do this. Specifically...

What exactly is a bin? Do the upper and lower boundaries of a bin correspond to the frequencies of two consecutive notes? In other words, because an 88-key piano has 7 octaves each with 12 unique notes, should I set n_bins=7*12=84 or equivalently bins_per_octave=7?  Or should several bins correspond to a single note interval?
Is fmin supposed to be the deepest note on the 88-key piano, ie the A note with a frequency of about 27.5 Hz?
Why do we need fmin? Is this some sort of reference point, similar to the equation from amp to decibels?
What are the difference between n_bins and bins_per_octave? Which one of them is better to use? I ask this because a research paper here (line 13 and 53) uses BOTH n_bins and bins_per_octave to analyze 88-key piano music. Why do they use both parameters instead of just one?


texts: system identification with delayed desired signal

texts:
I have a white noise data for 10 sec ( with 2 min silence or zero data at start and end)that is captured by sensor after delay t1( which I dont know). same data is captured after delay t1 and t2 by sensor 2. I want to compute transfer function(system identification) between sensor 1 and sensor 2.
when I use LMS( I know unknown is IIR) to adapt, I get filter which has zeros outside the unit circle, and impulse response which is non symmetric. I am told that this is typical of non-causal filter.

If I use sensor 1 data, capture portion of say 5 sec of white noise data from the first non zero sample, treat it as my input , and then do the same for sensor 2 data , treating it as desired signal for LMS, would this be a right approach?
Alternatively, assuming a group delay of filter_len/2, where filter_len is the length of LMS filter, sensor 2 data would be delayed by this amount , with respect to sensor 1. So, capturing white noise data starting from first non zero sample would amount to adapting the filter too soon( disregarding the group delay) and hence will give anticausal filter?
I tried capturing sensor 2 starting from group delay samples before the non zero sample, and repeat the LMS adaption, but the filler still looks all over the place.( see Figure attached , with LMS of filter length 256, step size 0.08)
Any suggestions on how can I solve this problem?
Thank You


texts: How can resonance be added to a fractional order Butterworth low pass

texts:
I am interested in how to add resonance (Q) to the magnitude response of a Butterworth low pass when it is expressed in the form:

$$
G^2(\omega)=\frac {1}{1+\left(\frac{\omega}{\omega_c}\right)^{2n}}
$$

I know how to add resonance $Q > 0$ for $n >= 2$ i.e. include it in one (and only one) of the second second order terms (if the equation is expanded as a product of second order sections):

$$
H(s) = \frac{ \omega_0^2 }{ s^2 + \frac{ \omega_0 }{Q} + \omega_0^2 }
$$

However, I wonder whether there is a general formula for resonance $Q$ that is a function of $\omega$ that can be used for fractional $n$ such that it can be added to the transfer function $G$ or $G^2$ to yield a new magnitude response?


texts: Feature Extraction from Accelerometers for classification

texts:
I am looking into algorithms/techniques for feature extraction for accelerometer data - to then be used in classification.

FFT is incredibly popular, however would this work on accelerometer data?  I've also seen PCA, but that isn't quite so popular.

Will PCA work as a possible feature extraction algorithm?


texts: Natural cubic spline interpolation versus cubic BSpline interpolation?

texts:
An answer here seems to shows the algorithm Mathematica uses to compute:
Interpolation[{{x1,y1},{x2,y2},{x3,y3},{x4,y4},{x5,y5}},Method->"Spline",InterpolationOrder->3].
I used to assume Mathematica found a natural cubic spline interpolation, but it seems Mathematica actually finds a linear combination of BSplineBasis functions that interpolate the data. What are the pros and cons of doing this using BSpline interpolation versus natural cubic spline interpolation? Also, what is special about the interpolation found using BSplineBasis functions?


texts: Educational purpose - What is the correct way to simulate a multipath fading channel which has ISI

texts:
I am following the codes for rayleigh fading channel given in the webpage: http://www.raymaps.com/index.php/m-qam-bit-error-rate-in-rayleigh-fading/

I am confused and need help in clarifying some basic questions regarding how to simulate a Rayleigh fading channel with $L$ taps? This is not clearly mentioned.

Question 1:  I have studied that a flat fading channel is one which has 1 tap. In the webpage, a flat fading is simulated by h=(1/sqrt(2))*(randn(N,1)+1i*randn(N,1)); where $N$ denotes the number of bits or the number of data points. So, there are $N$  taps and definitely this is not a 1 tap channel. Please correct me where wrong.

Question 2: Also, how to simulate a fading and flat fading with say $L = 5$ taps?

This is how  I have done but I am not sure how to properly do such that there is inter-symbol-interference as well. If number of paths, $L$ through which the signal reaches the receiver with some delay among the different paths. If the delay spread is greater than the symbol duration then we get inter symbol interference. In the following code, I am not sure if inter-symbol-interference is there or not and if I am correclty applying the theory to simulate the fading channel.

clear all
% Number of data points
N = 100;
L= 5; % number of taps
%Data symbols of 0/1
x = randn(N,1)>0;
%channel coefficients as Rayleigh random variable with variance 0.5 for real and imaginary component
h = 0.5*(randn(L,1)+1i*randn(L,1));
%Channel is an FIR filter
y = filter(h,1,x);


texts:
This task comes from tracking object on a steady background. So far I was able to remove the background and obtain binary masks like this:

I need to get bounding rectangle of the toy without the cord attached to it (on some frames it's quite thick). Sounds easy enough but I can't think of a good way to do it.

I tried cross-correlation with template of the toy I cropped. It wasn't really nice, plus I needed rectangle to be able to change it's size. And the cord added some mess.

meanShift() also didn't show good results because of the cord and occasional noice.

For now I stopped on findContours() with some tweaks, but it still produces some weird boundaries on what seems like okay image.

I keep thinking there must be some more simple and effective solution.


texts:
I was scanning the airwaves using my RTL-SDR and found a signal being broadcast locally that I was able to receive. When I configure the SDR Console V3 software to NFM, I get the following audio. (Wav file)

So I started investigating and here are my findings:

What I know:

According to my local FCC equivalent the license documents list it as coming from an electrical switch yard.

There is are two different frequencies for Tx and Rx.

The signal runs continuously.

What I assume:

That this is a 2-fsk/BFSK signal, due to the two tones it uses.

That it is broadcast in NFM. That mode does sound the best however, and lines up with licensing documents.

That it is actually data, not just some tone control system or random noise.

What I have done

In SDR Console V3 I created an IQ recording. (Large File)

I then created a flow graph in GNURadio to record a file using a sink file. (Large File)

I created a GNURadio flow graph using that file as an input to low pass filters and a FM Demod.

The Demodulated file was then opened in Inspectrum which produced the following image.

What I think I need to do next

I think I need to recover the clock of the signal somehow. I have had a look for a preamble but can't see anything that might match.

I thought it might be Manchester encoded but I can see long runs of 0's or 1's so I don't think it is.

I thought it might be clocked to the frequency of the local power supply given that it is in a switch yard, but it doesn't line up completely with rise or falls of the signal.

Am I on the right path with this?

Is there a way to determine the start of the message? The receiving radio has to be able to lock on to the start of the message somehow, right?


texts: Image processing for extremely small feature detection

texts:
I'm new to the world of image processing, and I need some hints to analyze an image.

I have a sensor which is made up of a high resolution CMOS camera sensor and there's a mask with a micrometric hole on it.

A light shines through the mask and the pinhole creates a circular projection on the camera's CMOS. Since the projection affects only 10-20 pixels out of the 8MP of the sensor, what would be a smart way to process the image the fastest way possible? I'm now using Matlab's native code but I want to code it myself.

The sensor is self-powered, so CPU processing is an issue.

Someone suggested wavelet analysis, I was wondering if that's really my best shot.

What do you think?


texts:
I experiment with the coefficient quantization in an IIR filter. I change some values like the passband ripple and also how many bits i want for the quantization method. In the first picture i have 1db passband ripple , 'round'and 7-bits as the quantization method:

In the second picture i have 1.3db as the passband ripple ,'round' and 6-bits as quantization method.

These are the PZP diagrams. I dont understand how i can extract useful information from the diagrams in order to say which design is better or to make comparisons between these 2 designs.How the Poles and Zeros affect my decision for what is better?

edit: the images becomes bigger and more visible when you copy the url and open or select view image. Also with the green is the quantized and with blue the un-quantized coefficients.

edit2:As Matt asks for magnitude frequency response I have add the extra information:

1db passband ripple , 'round'and 7-bits

1.3db as the passband ripple ,'round' and 6-bits

EDIT: After the answer from Matt i need to add more information:

I want to find a set of quantized coefficients that meet the specification with the minimum number of bits.
 My specification is an 8th-order IIR filter with the following transfer function magnitude response:
\begin{align*}
 0\,\mathrm{db} \pm 1.5 \,\mathrm{db} & \quad \text{for } \quad 0.2<|v|<0.3 \cr
 \lt −60\,\mathrm{db} & \quad\text{for } \quad |v|<0.14\text{ and } |v|>0.36
\end{align*}

I use this matlab function:  $ [b,a]=ellip(n,Rp,Rs,Wp) , $

after some trials i found that (the filter is implemented using second-order sections) :  $ [b,a]=ellip(4,1,60,2*[0.2,0.3]) $  and using 7-bit quantization i take the diagram above that i think verify the specification (if i make a mistake please correct me) ,

but if a use: $ [b,a]=ellip(4,1.3,60,2*[0.2,0.3]) $ and 6-bit quantization also i take a reasonable result. My concern is if my thought and procedure is correct or i have a mistake and i need to stick with 7-bits quantization?


texts:
Question: Can anyone refer me the Phong reflection model parameters for a face image taken for web-cam?

Details:

I am doing 3D reconstruction of 2D images using 3D Morphable Model as in this paper Vetter & Blanz, 2003.

There is a cost function which needs to be optimized using Newton Method. I need to give initialization values for the Newton method.

The below equation is the phong reflection model of an image. The ambient, diffuse and the specualr components are shown.

Phong Reflection Model wiki

Phong Parameters: Ks, Kd, Ka, α, Lm, V

I am estimating the above parameters using Newton method. The image is of human face and the camera is web-cam.

I gave the below values: Ks = Kd = Ka = 0.6
α = 5
Lm = [0 1 1]
V = [0 0 1]

Can anyone comment on the above values and what can be a good initialization values for the parameters??


texts:
I'm trying to understand how speech recognition works, when using HMM. I know this post is very long, but that's mainly because I'm a complete beginner, and I since I have no idea what is clear and what is not, I try to write everything as clearly as possible (which takes a lot of space). Help would be really appreciated.

I will be using O for observation sequence, X for state sequence, M for a HMM model with defined parameters.

My understanding - please correct me where I am wrong.

If I understand it correctly, the 3 basic operations with HMMs are:
1, Evaluation - how likely is observation O to happen, given a HMM model M - P(O|M). This is done by using the forward algorithm (dynamic programming).

2, Classification - I'm not sure I understand what this one is - is it finding out which state sequence is the most likely one to cause the observation O, given a model M? If yes, then it it is again done with dynamic programming, very similar to the forward algorithm, except you will be taking maximums instead of sums.

3, Training - You change the parameters of M, to raise the probability of an observation O happening with your given model. The forward-backward algorithm can be used for this, and it is proven that it finds a local maximum.

The actual speech recognition, here in the context of isolated word recognition:

For each word you want to recognize, you have a specific HMM model - here, the example I will use is that if you want to recognize 'one', 'two', 'three' you will have 3 HMM models M_1, M_2, M_3, each with say 3 states and 3 emmissions (that might be enough for 3 words?).
You make training samples of the specific words, and you train the models with these sounds - you do this, by  breaking up each utterance/sound of a certain word into small parts (~10ms), to each part you map a certain vector, with the vector being based on some spectral quality of this 10ms part - this is the observation at a certain time.

Now, if you want to recognize a word, you say the word, and you pick such a model M_i, that it is the maximum of P(O|M_i), i element of L, L being the number of words (and specific HMM models) that you have trained.

What I don't get:

What exactly can the vector/observation be? Can it be a relatively simple function, say approximation of the spectral envelope? Or is it usually a vector of R^n?

How is it, that if you say "ooone" instead of "one", P("ooone"|M_1) will still be the highest of the three trained HMM models? Is it because if the model was trained with 'one', the model already favors the observation that happens when the sound of 'o' is said, and the state transition of x_1 -> x_1 (where x_1 is likely to emit this observation) is also favored? Basically, is there some convincing proof/intuition, that pronouncing a certain word very slowly at certain parts will not cause it to be recognized incorrectly?

In the case of isolated word recognition, does the classification HMM operation happen at all? All I see is training, and evaluation (for finding the maximum of P(O|M_i).

Lastly, could some please explain how continuous speech recognition differs from this? Apparently, unlike the above, classification happens at some point - I don't really understand how it can be used. I can imagine that maybe for each phoneme, a certain HMM model will be used, and then in a long sentence, this sentence will consist of many of these models tied together.
How many HMM models will there be, how many states will there be? How will it be done, that a certain state sequence will define the words in the utterance?


texts: How can i measure shape elongation feature in matlab or python opencv?

texts:
How can i measure shape elongation faeture in matlab or python opencv?
There is a function library or i have to implement to my own?
In case there is no function, i found this definition:

"The less-common elongation shape factor is defined as the square root of the ratio of the two second moments in of the particle around its principal axes"

f_elong = sqrt (i2/i1)

Where i found these first and second moments?


texts:
In recent years I've started doing work in image analysis and have used the Barbara and Lena images. Early on I became aware of the interesting history surrounding the Lena image and that got me wondering: Who is Barbara, where was the image shot, and what was the original intent (presumably it wasn't created for image research)?

Many people cite Allen Gersho's lab or Marco Schmidt's (now defunct) test image database. However, this is a bit unsatisfying to me.


texts:
Let's say I have a decimal like 0.9

To convert to Q0.27 I just multiply by 2^(27)

But with Q5.27, with 1 bit representing sign, how would I use the full range of bits available to me, to say convert 0.9 to a number between 16 and -16?


texts:
Lets assume a pipeline X:

            in_buffer1 * window -> fft -> car_to_pol -> pol_to_car -> ifft -> out_buffer1
          /                                                                               \
in_signal                                                                                out_buffer1 + out_buffer2  -> out_signal
          \                                                                               /
            in_buffer2 * window -> fft -> car_to_pol -> pol_to_car -> ifft -> out_buffer2

where in_buffer1 and in_buffer2 overlap in time by exactly 50% and the window is i.e. a hamming window with an offset matching the buffer.

the out_signal should be just a reconstructed in_signal with a delay depending of the buffers lengths with no frequency leaks.

However lets then consider a following pipeline Y for a single buffer with two sifnal sources and one result

in_signalA(n) ->  in_bufferA(n) * window -> fftA(n) -> car_to_polA(n)  ->     imgA(n)     ->   pol_to_car -> ifft -> out_buffer(n)
                                                                       \                   /
                                                                        realB(n) * realA(n)
                                                                       /
in_signalB(n) ->  in_bufferB(n) * window -> fftB(n) -> car_to_polB(n)

Then apply the same overlapping oerlapping as per the pipeline X where the n is the number of the pipeline, i.e. 0 - no offset, 1- 50% offset,
This should produce a 'vocoder' effect.

However I do get cracks on the output signals which suggest there is a frequency leak I'm not compensating for. Any idea what I'm missing here?


texts: Why Is the Total Time Equal to $ N \cdot {T}_{s} $ and Not $ \left( N - 1 \right) \cdot {T}_{s} $ In the Context of DFT?

texts:
In the definitions of the DFT
DFT
$$
X(j)=\sum_{k=0}^{N-1} x(k) \exp \left(-i 2 \pi\left(\frac{j}{N}\right) k\right)
$$
Let us say, if we have $10$ points, $N=10$, each sampled at $0.2$ seconds, why is the total time used for calculating the frequency resolution is equal to
$$
\frac{1}{N\Delta t }
$$
where $k$ will run from $0, 1, 2, \ldots , 9$.
If the first point was at zero, the sampled time will be at $0.2$, and the last sampled point will be at
$$0.2\cdot (N-1)= 0.2\cdot 9 = \mathbf{1.8 \ \rm s}$$
Rather the total time is equal to $0.2\cdot N= 0.2\cdot 10=\mathbf{2.0 \ \rm s}$ in the frequency step.
P.S.: I have seen the query and the discussion How do I measure the time duration?

Here $\Delta t = 0.1 \ \mathrm{s}, N= 11 (\text{Eleven data points}),
k= N-1$; So
$$\text{total signal duration} = k\cdot \Delta t = (11-1)\cdot \Delta t= 1 \ \text{second}$$
This agrees with the $10 \rm \ Hz$ sampling rate, i.e., $10$ points were collected in $1 \ \text{second}$ and the $11^{th}$ point belonged to the next cycle.


texts: Possibilities to get the amplitude of a signal at specific frequency from IQ

texts:
I'm newbie to the dsp topic and my understanding of math symbols is limited.
So at the moment I'm little bit lost to find the right resources.

I receive an IQ stream from a SDR.
Sample rate: 250k/s
The signal I'm interested in is at -10kHz.
I need the amplitude of only the -10kHz frequency.

I've found out how to shift the IQ data frequency up/down. But I thought -10kHz is a nice place to be away from DC, and for buffer sizes and sample rate.

What is the most simple solution? Are there other possibilities next to one of the Fourier transformations? Looks like overkill for my simple signal.

But I've found lot of theory stuff handling IQ signals with only one frequency, but nothing enlightening to understand the demodulation of one frequency out of a mix of many from IQ.

I know about the phase change. But something I miss to understand how to collect and calc the right data for my frequency from the IQ.

What would be a nice book (German/English) which explains the stuff without getting lost in math vocab and theory after a few pages? Anything from a software developers perspective?

Thx.


texts: Study materials for Decision feedback equalizer

texts:
Looking for books, course notes/slides, or some tutorials on decision feedback equalizers with in-depth details.
I appreciate your help.