texts: A Laymans Fourier transform and harmonics explanation?

texts:
Please bear with me, I know some of you will scoff but I have looked on Wiki and in various literature (see below) and can't quite get a handle on a few things.

I am a general business programmer, in so much as I very rarely perform any complicated maths/physics related programming.

Anyway, whilst at Uni and studying Computer Networks I was introduced to Fourier analysis down at the physical level. I had trouble understanding the actual application of the analysis (as described in A Tannenbaum's  Computer Networks) he just throws in some equations and asserts that they are used to approximate a signal modulated to carry data.  I could also never understand the relationship between harmonics and the signal.

Lastly, in summary, is this analysis used to pick up from a physical medium, ie re-construct, retrieve the data initially sent.

Thanks,


texts: necessity of patch size in edge detection algorithms

texts:
I have recently started learning about the image edge detection methods. I studied few of these methods including Sobel and Canny. I observed that edge detection is performed using patch/kernel size as 3x3.

My question is, Is it necessary to fix the patch size as 3x3? OR can we take it as 5x5, 7x7, 9x9 etc.

I have already studied few technical papers that says that patch size is taken as low as possible in order to make computational complexity as low as possible


texts:
I have a camera (iPhone), I have a 3D control object in the image that I know its properties very well. (My control object). There is also a secondary object in motion. The ultimate goal is to establish the 3D trajectory of the moving object for a given period of time. (Tracking)

I like to ask, could I find out?

Phone's distance to the control object (for sake of discussion, let's assume camera is at certain height and certain distance neither of them is known but camera is perpendicular to the surface which is known)
The secondary object where I can locate the object in each subsequent frame.  My goal is to estimate its 3D trajectory as I have indicated above.

Bonus question, we can make the system such that the phones distance to the control object can be set (although not preferred), would this help me with the second point?


texts: How do I generate a QAM signal using Python?

texts:
I have created a 16-QAM Constellation diagram ranging from {-3-3j ... +1+1j}. I then transform a given input string into symbols where each symbol is log2(16) = 4 bits wide. I then map each symbol to a constellation point, i.e. symbol 15 would be the constellation point 1+1j and symbol 0 would be -3-3j
The next step is to create a playable sound, for this I am trying to generate samples using the following formula that I found on Wikipedia and in other places:
s = sin(2*pi*fc*t)*I(t) + cos(2+pi*fc*t)*Q(t)

sample = np.math.sin(w * t) * constellation[symbol].real + np.math.cos(w * t) * constellation[symbol].imag

Where w = 2 * pi * fc
And finally, appending the single samples to a full sample list and multiplying by the maximum amplitude (32767):
samples.append(int(sample * 32767))

before saving the data into a .wav file:
qam_file.writeframes(struct.pack('%dh' % len(samples), *samples))

Something doesn't add up, and I have no idea what, am I using the correct formula or am I doing something wrong with the In-phase and Quadrature Components?
Any help would be greatly appreciated since I am fairly new to DSP.


texts: What happens to crosscorrelation of two noise generator if one of them were to be inverted

texts: I have two noise generator u1 and u2 and their crosscorrelation is given, what will happen to this crosscorrelation if u2 or u1 were to be inverted?


texts: Use MATLAB to Restore a Signal from a Given Degraded Signal Using Tikhonov Regularization

texts:
Anyone could share how to develop an application in MATLAB to restore the signal from a given degraded signal using Tikhonov regularization i.e restoring the signal $f$ via solving

$$
\min  || g - f \circ h ||^2  + \lambda ||\sigma \circ f ||^2
$$

where $\sigma = [1,-1]$. and $\circ$ denotes the circulant convolution operator. The setting of the value epsilon is dependent on the signal to noise ratio.

Given an audio file with noisy convoluted measurement,

The following routine are expected
Tikdeconv.m

Tikdeconv De-convolue a signal using Tikhonov regularisation J=deblurw(I,PSF,LAMBDA) de-convolute input signal $I$ by PSF using Tikhonov regularisation with regularisation parameter LAMBDA returning de-convoluted signal J. The assumption is that the signal $I$ was created by convolving a signal $J$ with kernel PSF and possible by adding /guassian white noise. You may not call or use any code from the following built-in function in MATLAB:deconvwr,deconv,deconv2

Any feedback on how to solve this would be appreciated. Thank you


texts: Wavelets: Reconstruction Filters for 2 Level Decomposition

texts:
I am trying to understand the reconstruction part after discrete wavelet decomposition and how do we get approximations and details at various levels. Most textbooks show the complete reconstruction diagrams but not how to get an approximations and details at various levels.

Let us say we started with a s[n] and DWT is performed at level 2. I am following the diagram in Fugal's Conceptual Wavelets.

We pass a high pass filter, we get (cD1) and a low pass filter, we get (cA1) on the signal. The cA1 is downsampled by 2, and then passed through the same high-pass and low pass filter giving us detail coefficients at level 2 (cD2), and approximation coefficients (cA2).

How do get approximations A1, A2, and details D1 and D2?

a) If we upsample cD1 by 2 and pass a reconstruction filter on it we will get D1.

b) If we upsample cA2 by 2, and pass a low pass reconstruction filter, we should get A2.

c) If we upsample cD1 by 2, and pass a high pass reconstruction filter, we should get D2.

d) How is A1 obtained? As per the diagram in Fugal (Conceptual Wavelets), he combines the output at level 2, upsamples it again, passes it through a low pass reconstruction filter and gets A1. I believe this could have been done by passing the low pass reconstruction filter directly on cA1 as in point (a-c).

I am asking because MATLAB displays the approximations and details at various levels.

Thanks.


texts: pdf of the sum of gaussian distributions using fft Ask

texts:
I am trying to derive the pdf of the sum of independent random variables. At first i would like to do this for a simple case: sum of gaussian random variables. I was surprised to see that i don't get a gaussian density function when i sum an even number of gaussian random variables. I actually get:

which looks like two halfs of a gaussian distribution. On the other hand, when i sum an odd number of gaussian distributions i get the right distribution:

below the code i used to produce the results above:

import numpy as np
from scipy.stats import norm
from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
%matplotlib inline

a=10**(-15)
end=norm(0,1).ppf(a)
sample=np.linspace(end,-end,1000)
pdf=norm(0,1).pdf(sample)
plt.subplot(211)
plt.plot(np.real(ifft(fft(pdf)**2)))
plt.subplot(212)
plt.plot(np.real(ifft(fft(pdf)**3)))

Could someone help me understand why i get odd results for even sums of gaussians distribution? Thanks in advance


texts: What are features in signal processing using a more tangible or less abstract description

texts:
I'm learning about signal processing for music, and in my reading I keep finding the term "features". I find that they never really give a concrete example of what they mean by "features". They'll say stuff like features related to timbre, but they won't really go into further depth than that. I also took an image processing class once and found that I didn't really understand it there either. I understand what some literal features of a sound would be, like say distortion, but I have a feeling that this isn't really what they're talking about when they say features.

Example
I'm reading something on Dynamic Time Warping(an algorithm that seems similar to local/global alignment) and they say this

To compare two different features x, y ∈ F, one needs a local cost
measure, sometimes also referred to as a local distance measure, which
is defined to be a function c :F×F → R.

My understanding of Dynamic Time Warping is that two similar waves(∿), aligned to each other by measuring the distance between each subsection of one wave to each subsection of the other wave and then choosing the alignment that minimizes the total distance between all the wave sections. I'm doing this for soundwaves, so why are they describing it as comparing two different "features"?

Obviously it can't be described in a way that is 100% tangible, but I'm looking for a description of features that one would be able to draw out or make a picture of.


texts: $n\textrm{ Hz}$ waveform sampled at $m\textrm{ Hz}$ per second

texts:
Here is an example of plotting a square wave given in SciPy Documentation

A $5\textrm{ Hz}$ waveform sampled at $500\textrm{ Hz}$ for 1 second:

from scipy import signal
import matplotlib.pyplot as plt
t = np.linspace(0, 1, 500, endpoint=False)
plt.plot(t, signal.square(2 * np.pi * 5 * t))
plt.ylim(-2, 2)**strong text**

The documentation page states that

  The square wave has a period 2*pi, has value +1 from 0 to 2*pi*duty
  and -1 from 2*pi*duty to 2*pi. duty must be in the interval [0,1].

Why in the example 5*t is multipled by 2*np.pi?
Doesn't the signal.sqrt() function
take care of multiplying the frequency 5*t by 2*np.pi?
What is the meaning of sampled at $500\textrm{ Hz}$ for second? I understand it is a $5\textrm{ Hz}$ waveform since 2*np.pi*5*t takes care of creating $5\textrm{ Hz}$ wave. But I do not understand

  sampled at 500 Hz for 1 second


texts: How to determine the step response given a transfer function?

texts:
I'm new to Signal Processing. I would like to know how to draw the step response given a transfer function.

For example, given $G(s)=\frac{(1-\frac{s}{3})}{{(\frac{s}{0.1}+1)(\frac{s}{100}+1)^2}}\cdot 10$ a transfer function, I calculated the initial value with the "Initial Value Theorem": $g(0^+)=\lim_{s\rightarrow \infty}sG(s)=0$, but other than that, I would also like to calculate the final value and the settling time.

I've been also given the Bode diagrams (I don't know if it helps):

I think $G(s)$ is a transfer function of the third order, so do the second-order laws also apply there, where the settling time is calculated as $t_{s}=\frac{4}{\zeta \omega_{n}}$?


texts: Missing known sinusoid in FFT of sensor data


texts: Why OFDM is more appropriate with MIMO than other multiple access method?

texts:
As read from MIMO-OFDM wikipidia page says

  Research conducted during the mid-1990s showed that while MIMO can be used with other popular air interfaces such as time-division multiple access (TDMA) and code-division multiple access (CDMA), the combination of MIMO and OFDM is most practical at higher data rates.

Why does OFDM has more higher data rates and they combined more advantaged? Thank you


texts:
I have a signal with a sample time of 0.5 micro seconds and I would like to shift this signal by a fraction of the sample time, say by 3 nano seconds.

I have read a few online resources about fractional delay filtering and about using the FFT and IFFT to perform such delay. Can someone point me to some theory on this or give me an idea on how to implement it.

For the regular shifting of the signal for integer samples, I have implemented this by shifting the signal by the required number of samples and adding zeros at the beginning. Is this approach correct?


texts:
The transform function of a digital FIR filter can be characterized by
$$
H(z) = \sum_{n=0}^{2N} h(n)z^{-n}\quad\text{where $2N$ is the order of the filter.}
$$

But when we use $2N$ or $N-1$ , $2N+1, \ldots$, on any basis chosen, is it related to linear phase or  filter characters , or related to various configurations (low pass, high pass, band pass and band cut).
On the other hand , is $h(n)$ the coefficients or weights ?
Can you give me a realistic example of bringing the idea closer?
Are they converted to electrical signals later? Or what?

I Need Clarification Please.


texts: Adaptive filter to scale and phase shift two sensors output

texts:
One way of separating downgoing and upgoing wavefields in offshore seimic processing is to add signals from hydrophone and vertical component of the geophone (they are co-located). Hydrophone only registers a change in the pressure whereas geohphone as well as registering a change in seismic field also reacts to the direction of arriving wave, ie. if it is coming from above or below (sensors are deployed on sea floor). But the scale of the two are different as they record two different physical quantities. In addition, there is a phase difference between the two. To do the summation, I first need to re-scale one to the other and phase shift. The phase shift is slightly frequency dependent. I would like to make a filter that does the above on a period of data and then use the filter on all incoming data. The figure below shows an example of an upcoming signal (hydrophone signal has been downscaled for the exmple. In reality it is of much higher amplitude). If the two were of same scale and in phase then summing them will double the signal's amplitude while reducing the noise coming from above. I have some solution in time domain (assuming constant phase shift and scaling factor) and frequency domain by multiplying one's spectra with the transfer function between the two (found from spectral division). I have also used 'nlms' adaptive filtering. But some noise still remains, enough to make trouble later on.
Wonder if anyone has a suggestion on how to make the filter or some Matlab code.
Thanks


texts: Viola Jones algorithm: Haar like features, how are the features extracted?

texts:
If I have an image like this
1 2 3 4 5 6 7 8
a b c d e f g h
.....

And I apply a Haar like feature with template
1 1 1 1
-1 -1 -1 -1

Then in the first position we get X1 =  1+2+3+4+a+b+c+d
If we slide one side to the right we again get X2 = 2+3+4+5+b+c+d+e.
This way we will get X1 and X2 and X3 and so on. Now how are these values combined to get the feature? Because when we say a feature we are not just running that template in one place, rather we will run it over multiple places in the image. It gives lots of values like X1,X2 and X3 and so on...now how are those combined to get the final feature which will be passed to Adaboost?


texts: Trying to implement matlab pwelch function in python using scipy welch


texts: Finding the frequency response of a filter defined as a Z-domain transfer function

texts:
I am software engineer studying DSP on my own. I came across this question in a graduate textbook and am not sure how to proceed.

Determine the frequency response for the filter:
H(z) = 1.5 / (1 + 3z^-1)

Clearly, this is a z-domain transfer function. How do I get the frequency response? I am thinking of converting it into a time-domain difference equation, calculating about 128 points, running it through MATLAB's FFT, plot it, and voila, the final curve would be the answer. Is that a right approach? Is there a more analytic/math-y answer?


texts: Two generated tones creates a side tone - how to avoid the side tone

texts:
I'm trying to create two tones at different frequencies. So I have one tone at 15000 Hz and one at 15400 Hz using Audacity. The problem i'm running into is there is a third tone created at a much lower frequency I can tell exists from both my ears and a spectrogram.

How can I get rid of the lower tone? Where does it come from? I only have two tracks


texts: Problems with IFFT not being symmetrical

texts:
I have two signals, a measurement and a reference which I have performed an FFT on.  They have both been windowed with a Hanning window, and now I would like to deconvolve these to get the impulse response.

I have followed these steps:

Divide my reference FFT by my measurement FFT. (Point division of complex numbers).
I then take the conjugate of this resultant array.
I then perform a forward FFT on this.
I then take the conjugate again.
I then divide by the length to normalise the signal.

If I am to understand correctly, this should give me the symmetrical deconvolution of the two signals.  However, for my input of 1024 samples, I get the first 512 samples equal to 0, and then the second half of the samples are a little bit noisy with a definite impulse at the expected delay time.

My measurements are actually completely electronic so should be noise free.  I can confirm that my two input FFTs are symmetrical (I am using the full FFT, not just the positive values), I am not doing any 'FFT Shift'.

If anybody has any suggestions or hints, I would greatly appreciate it.

Thanks


texts: How significant is the time lag between inertial IMU data in a inertial analisys?

texts:
I'm working on a MCU system that use an IMU (ICM-20948) which has two different chip inside: one provides accelerometer and gyroscope data at $1125$ Hz (configurable at $125$ Hz), the other one magnetometer data at fixed $100$ Hz (not configurable).
Now I'm wondering how much the non-synchronicity of the data could be a problem in a motion processing algorithms.


texts: Accelerometer Bias estimation with kalman filter

texts:
I have to estimate biases of a 3-axes accelerometer by modyfying an existent kalman filter mounted on a drone. The biases are assumed constant.
The filter has 9 states: position (xyz), velocity (xyz) and acceleration (xyz) in navigation frame and no inputs.
Measurements are taken from GPS, IMU and baro.
An ahrs filter takes gyro, Imu e mag measurements to estimate roll,pitch and yaw.
To include Bias in the state Vector, I considered the measurement of accelerometer as :
Acc = acc_ned + Cbn * Bias;
The rotation Matrix from navigation to body frame Is needed since my Bias Is in body frame.
With this equation I wrote the relative modified matrices A and C and applied the Linear discrete-time kalman filter. The estimated Bias is substracted in feedback to the accelerometer measurement entering the ahrs filter.
I can't maKe this loop to work; my biases don't converge to constant values.
I'm trying simulations on MATLAB/simulink.
Has anyone ever tried to estimate accelerometer's biases with kalman filter? Has anyone any ideas?
Thank you in adavnce


texts: Convolution of delta function in discrete time with parameter

texts:
If I have a delta function of the following type

$$ n\delta(5n-20) * 3^{5n}u[n] $$ How can I calculate this convolution?

I have thought of the property $$ \delta[n-n_0] *x[n] = x[n-n_0] $$

Then if I say $$ k=5n $$ and $$ t[n]= 3^{5n}u[n]  $$

It should be that the convolution is equal to :

$$ nt[\frac{n}{5}-20] $$ and replacing t $$ n 3^{5(\frac{n}{5}-20)}u[\frac{n}{5}-20]$$

Is this right? Or am I wrong somewhere?


texts: Equalizing audio under linux / Raspberry Pi

texts:
I have an audio streaming receiver written in Python running on a raspberry pi. This receiver writes the streamed audio out to an alsa device, which in turn passes it to an amplifier which is driving a DML speaker that I made. DML speakers require rather heavy equalization to make them sound correct - so I want to add in digital signal processing on the rpi.
Using Room EQ Wizard, I have created a correction filter. I am able to see the raw equalizer values (frequency, gain and Q), and also export the filter as an impulse response wav. I have tested this IR wav using shairport-sync (I even wrote the directions) and it works very well.
The question is how to implement the equalization under linux (headless, no GUI) for audio that is coming from a python script.
I think there are three different routes that could work -

using a software based graphic equalizer
applying convolution to the audio signal either under ALSA after the audio is written by python
or directly in python before it is written

There is an alsa equalizer plugin, however this is too limited as there aren't enough bands, the bands are preset frequencies, and the Q is not adjustable.
There is also an experimental alsa convolution plugin, which unfortunately I can't get working.
I have also found this question which lead me to ponder whether or not python could apply the convolution directly. Unfortunately I am not a python (or audio processing) expert, so this would take quite a lot of effort to figure out.
Is anyone aware of an equalizer or alsa plugin that will do what I need?


texts: How should I modulate a $16$ bit input using FSK?

texts:
I need to graph a $16$ bits binary sequence, assuming that it'll be a FM analog signal I need to apply FSK to it.

So, I think I should use the formula

$A \cos(2\pi f_1 t)$ for binary $1$'s

$A \cos(2\pi f_2 t)$ for binary $0$'s

But I'm really confused about how to set a method for getting the graph for my $16$ bit sequence ...


texts:
Question:

Is there an algorithm that does (what I would call) "localized" gamma correction?

That is, applying different gamma corrections on different parts of an image, using a masking filter or something along those lines.

I have two images that seem to be correlated by a mix of exposure and gamma corrections, but not equally on the whole image:

First image

Second image

For further details on how these images were created, please refer to this (unanswered) earlier post.


texts: How to get histogram data object from matlab

texts:
Lets say I have a matrix x=[ 1 2 1 2 1 2 1 2 3 4 5 ].
To look at its histogram, I can do h=hist(x).

Now, h with retrieve a matrix consisting only the number of occurrences and does not store the original value to which it occurred.

What I want is something like a function which takes a value from x and returns number of occurrences of it. Having said that, what one thing histeq does should we admire is, it automatically scales nearest values according!

How should solve this issue? How exactly people do it?

My reason of interest is in images:

Lets say I have an image. I want to find all number of occurrences of a chrominance value of image.


texts: What is the $\mathcal{Z}$-transform of a constant?

texts: The Fourier transform of a constant exists. Can anyone please tell me what the $\mathcal{Z}$-transform of a constant is? Thanks in advance.


texts: Help me out in finding the problem in my given MATLAB code of Rectangular Contour Algorithm(RCA)

texts:
I am implementing the Rectangular Contour Algorithm (RCA) used for blind channel equalization on 32 Rectangular QAM Signal, polluted with noise and channel effects. It is described in:

  Khan Q. - Two Novel Blind Equalization Algorithms for Rectangular
  Quadrature Amplitude Modulation
  Constellations

I am not achieving reduction in the ISI value in dB as the iterations increases. When I use exact value of $\mu$ as given in research paper it diverges giving equalizer tap weights to NaN and everything goes to NaN. But decreasing the step size $\mu$ it becomes stable, whilst ISI increases instead of decreases.

Here is my MATLAB Code:

clc;
clear all;
close all
m = 32;
SNR = 30;
ITERATIONS=5000;
RUNS= 3;
%h=[-0.005-j*0.004 0.009+j*0.030 -0.024-j*0.104 0.854+j*0.520 -0.218+j*0.273 0.049-j*0.074 -0.016+j*0.020];

h=randn(1,7)+i*randn(1,7);
%h=h/norm(h);

N =length(h);                                      %equalizer length
M =32;
x = 0:M-1;
y = modulate(modem.qammod(M),x);
%32 rect qam
 A1 = [-7 7 -7 7 real(y(5:28)) 7 -7 7 -7]';
 B1 = [-3 -3 3 3 imag(y(5:28)) -1 -1 1 1]';
Y1=A1+j*B1;
 %scatterplot(Y1)
QAM = Y1';
Es = mean(abs(QAM).^2);                     % Average Signal energy
S=floor((N-1)/2);

nm = N-1;
Res_ISI = zeros(1,ITERATIONS);
MSEn = zeros(1,ITERATIONS);
bn=zeros(1,ITERATIONS);
%h = h./sqrt(sum(abs(h).^2));
R = length(h)-1;
cl = length(h)-1;
varn = Es*10^(-SNR/10);
stdn=sqrt(varn);
stdn_by_2 = stdn/sqrt(2);

for runs=1:RUNS
    R_S=[];
    E_O=[];
    W=zeros(N,1);
    W(S,1)=1;
    C = zeros(N,1);
    X = zeros(cl+1,1);
    disp(['Run = ' num2str(runs)]);
    for i = 1:ITERATIONS
        x = QAM(ceil(rand*m));
        X = [x; X(1:cl,1)];
        Ns = (randn+j*randn)*stdn_by_2;
        c = h*X + Ns;
        C = [c ; C(1:nm,1)];
        YY = transpose(W)*C;%Equalizer Output
        E_O=[YY E_O];
        R_S=[c R_S];

%RECT Contour Implementation
mu = 1e-15;
A=7;
B=3;
R_rect=1.7927;%32 RQAM
P=1;
YR=real(YY);
YI=imag(YY);
T1=B*YR+A*YI;
T2=B*YR-A*YI;
S1=sign(T1);
S2=sign(T2);
S3=A*YR-A*YI;

%Del Jrect Equation
del_Jrect=((abs(T1)+abs(T2)).^P-(abs(A*B)*R_rect).^P).*((abs(T1)+abs(T2)).^(P-1)).*((B*S1+B*S3)+i*(A*S1-A*S2));
del_Jrect=del_Jrect.*conj(C);
%Taps Update Equation
W=W-(mu*del_Jrect);

TA=abs((h.*conj(W)')).^2;
TB=max(TA);
res_isi(1,i)=(sum(TA)-TB)/TB;
    end
    Res_ISI = Res_ISI+res_isi;
%    MSEn = MSEn+mse;
 scatterplot(R_S);
 %scatterplot(filter(W,1,R_S));
end
% end

Res_ISI=Res_ISI/runs;

%____________________________________________________________
 figure;
 plot(10*log10(abs(smooth(Res_ISI,100))));grid on;
ylabel('ISI(dB)'); xlabel('Iterations');


texts: whether the system is linear or not for the given problem

texts:
Given the system:
$$y(t)=x(t+1)+x(t−1)$$
is the system linear?

For a system to be a linear first it should satisfy zero input and zero output. How can we calculate output at 0 input if the system depends on future or past or both? Please explain with steps.


texts: determine background noise level via cross correlation

texts:
The setup:
Sound is being played by the speaker on the iphone. At the same moment i record sound via the microphone. Obviously the speaker sound will be recorded in addition to all the external noise.

My goal:
To determine the level of all external sound which is recorded in addition to the sound from the speakers.

(no implementation yet)
So far my approach is to cross correlate the input stream with the output stream to get the delay between both of them and an integer representing the level of similarity at this point, which will be used to estimate external noise. (if the int drops below a given value the noise is too loud)

Has anyone experience with stuff like that, if this approach will be accurate enough to make a guess if the external noise reaches a certain level of disturbance?
Or will i have t do further analysation in the frequency domain via fft for example.

I'm completely new to this kind of problem so be kind with me, if I'm totally wrong ;)


texts: Interpolating and Decimating CIC filter based on sinc function

texts: How to realize a sinc function using CIC filter for decimation and interpolation ??? Can I combine interpolation and decimaton methods inorder to get the complete response of a sinc function???


texts: Calculating Savitsky-Golay filter coefficients - normalization

texts:
Using the matrix equations found on Wikipedia I wrote some Python code to compute the Savitsky-Golay coefficients for $m$ data-points, using a polynomial of degree $k$. Here it is:

from numpy import *
from numpy.linalg import *

k = 4       # Degree of polynomial to use for regression
m = 25      # Window size, that is, the number of points. Must be an odd number.

z = range(-(m-1)/2, (m-1)/2+1)  # z is the distance from the central point.

jacobian, b = [], []
for i in z:
  for j in range(k+1):
    b.append(i**j)
  jacobian.append(b)
  b=[]

J = matrix(jacobian)

C = (J.T*J).I*J.T

I do realize that this code is in no way optimized and it could've been done in a much cleaner way, however I'm not a (professional) programmer nor a mathematician and would be really happy if I could get it to work at all.

My problem is that the coefficients computed become rather small, using the parameters above ($k=4$, $m=25$) on the order of $\pm 10^{-6}$. If I "go one step back" and multiply with the determinant I get really huge numbers ($\pm 10^{22}$).

I tried finding the greatest common divisor of these but got the result 2048 which leaves me with huge numbers (and it can't be right, result of numerical errors due to small/huge numbers or something?). I would like to have the coefficients in the form they're usually found in tables, with a normalization constant for each derivative, and an integer coefficient for each data point. Since I'm interested in the point of maximum curvature in a data series I want to calculate its 3rd derivative, I've been unable to find a table with these coefficients.


texts: Multi-channel filter for correlated noise suppression

texts:
Short question
We need to cancel the correlated noise from recorded datasets, as shown in the figure below, coming from arrays of charge sensitive sensors.
I would like to ask to this community to point me to the relevant literature about correlated noise suppression. In addition, I would like to know if there exists some code and/or algorithms, and a base textbook that can serve as reference in the future for offline signal processing (I believe it is referred as DSP).

Example showing 36 out of 64 total channels of a triggered event where charge signal is present. The upper pad shows the single channels where the amplitude of the waveforms are shown represented by the color scale. This view shows evidence for a correlated noise across chanels as vertical runninbg "valleys" and "cliffs". The lower pad shows the waveforms of the same channels overlapped and plotted with a different color each. The strong pulses correspond to electrons "collected" at our anode electrode (pixels) and are our "signal".
The longer story
I am analysing digitised data from several charge sensitive amplifiers (here a link) of a particle detector, consisting in a liquid argon time projection chamber (LArTPC).
As shown, the recorded waveforms are strongly affected by correlated noise/disturbances. I believe that this noise is due to some ground loop that we will identify and remove in some future setups of our detector (still in R&D stage) and to some EMI of nearby power electronic (and a train station closeby our laboratory).
In order to get some accuracy from the analysis of the charge (taken from the hight of pulses), we need to figure out how to filter out the correlated part of the noise among the channels. At the moment, we are subtracting from each channel either the mean or the median waveform calculated using all the channels. Both methods work not too bad when the SNR of the pulses is high as in the shown example. However, when the SNR is lower the performance really drops making the analysis almost impossible. Being able to study also small signal events is fundamental for us, as we need to investigate the evolution of the liquid argon purity achieved over time.
What I would like to develop is a filter based on the waveforms in the "No signal region" (see the figure), where no signal can fall, and apply it to the middle region where charge signals can occur.
I have digged quite a lot into the literature, but it was almost a random walk as I am not expert of signal processing/conditioning. Thus, I am not sure that I could find (or understand) what I was looking for. The only things I could find (and understand) are the following 3 papers: Trojanoski et al., Wang et al., Ozbek.
It seems that the first two references propose to estimate the Green function (in frequency domain) between two channels (R$_{ij}$) averaging noise only (seismic) waveforms over several acquisitions. This has the advantage of cancelling the uncorrelated noise for getting a good estimation of the correlated one. If I am right, this is not really applicable to our case as the ground loop and EMI induced noise is not really stationary across several events. The third link proposes an adaptive filter (in time domain) and assumes non stationary correlated noise. However, the credentials of my institution don't allow me to get that paper and thus I do not know whether it is relevant for our needs.
At the moment, I am mainly looking for filters that minimize the RMS (Wiener etc.) as I can better understand the math and what they do. It seems that adaptive filters also are interesting especially if the noise is not stationary. However, once we have managed to implement in our analysis a better filter, it would be good to make a next step trying to develop and implement a filter with minimum or very controlled signal distortion.


texts: The impulse response, $h(t)$, of a system

texts:
I have this differential equation
  $$
    y''(t) + y'(t) + y(t) = x(t) - 3x'(t)
    $$

I found the system response of the system:

$$
        H(f) = \frac{1 - 6\pi fj}{(2\pi fj)^2 +2\pi fj+1}
    $$
which is equivalent to

$$
        H(s) = \frac{1 - 3s}{s^2 +s+1}
$$

My question is: How to find the impulse response? If I factorize the denominator, it's not gonna be pretty. e.g., $s + \frac{1 +\sqrt{-3}}{2}$

Also, what does the following mean?

Find the output for $t \ge 0$ with initial conditions $y(0) = 1$, $y'(0) = -1$ when the input is $2cos(6t)u(t)$. Identify the zero state, the zero input, the transient, and the steady state responses.


texts: Why adaptive filter does not work in my application

texts:
I got a problem when I was trying to denoise a signal. Actually, it is a simple simulation. The signal is the addition of a step signal (The info I wish to get) and a sinusoidal one (the noise I wish to remove). See below
However, I tried different parameters of using the adaptive filter, it simply cannot filter out the sinusoidal noise from the step signal. See figure like this.

Any suggestions will be greatly appreciated!

Below is the matlab code

clear all
close all

%% walking induced noise
t = [1:5000]*1e-2;
f = 0.1;
WalkNoise = 1*sin(2*pi*f.*t);%+1.5*cos(3*pi*f.*t);
WN = WalkNoise + 0.05*randn(size(WalkNoise));

figure
subplot(3,1,1)
plot(t,WN);
title('Noise');

%% signal
h1 = 14; % height of the signal 1
h2 = 18; % height of the signal 2
L = 5000;  % total length of the signal
bp =2500;  % location of break point
x1 = h1*ones(1,bp);
x2 = h2*ones(1,L-bp);
Sig = [x1,x2];
Sig = Sig + 0.1*randn(size(Sig));

subplot(3,1,2)
plot(t,Sig);
title('Signal');

%% walking-induced-noise + signal
NoisySig = Sig + WN;
subplot(3,1,3)
plot(t,NoisySig); hold on
title('Signal + noise');

%% adaptive filtering
figure
plot(t,NoisySig); hold on
title('Signal + noise');

mu = 0.001;            % LMS step size.
ha = adaptfilt.lms(20,mu);
[y,e] = filter(ha,WN,NoisySig);
plot(t,e,'r');
legend('Signal+noise','Filtered using Adaptive filter');


texts: How to properly time-shift a function in relation to another function

texts:
I'm stuck on a time-variance question in my  homework and I think it has something to do with my understanding of how time shifting functions works.  The problem asks to determine whether the system defined below is time-variant or time-invariant.
$$ S: x(t) \to y(t) = x(t-2) + x(2-t) \tag{1}$$
I solved it such that
$$ x_{delayed}(t) = x(t - \tau)\tag{2}$$
and
$$ y_{delayed}(t) = x(t-\tau-2) +x(2-t +\tau)  \tag{3}$$

Then I delayed the original output $ y(t) $ from $(1)$ to
$$y(t-\tau) = x(t-\tau-2) + x(2-t+\tau) \tag{4}$$
which I believe is incorrect.  I think that instead, it should be
$$y(t-\tau) = x(t-\tau-2) + x(2-t-\tau) \tag{5}$$

I believe my problem comes from my understanding of how to delay a function related to another function.  Searching around didn't help much.  Ultimately my question is, should
$$ y(t) = x(-t) \tag{6}$$
$$ y(t - \tau) = x(-t + \tau\tag{7})$$
or
$$ y(t-\tau) = x(-t-\tau)\tag{8}$$
Using $(7)$ gives my my current answer, $(4)$, while using $(8)$ gives me what I believe is the correct answer, $(5)$.  I have found proofs for both answers which only confuses me more, so which is correct and how should I time shift the function? Any help is welcome and appreciated, thanks.


texts: Segmentation of high noise, high density and low resolution image

texts:
I've come across a particularly difficult segmentation task.
It's optical coherence tomography data, that has a lot of noise and, compared to the particle size (radius= ~3-6px), low resolution.
This, combined with a many overlapping particles, has made it impossible for me to segment out the particles with decent accuracy.
I've played around with thresholding, watershedding, morphological operators and h-dome transformation but couldn't come up with results that I've found satisfying.

The region of interest is already cropped out (inside of blue boundary):

This is what 10 frames of the recording look like:

Maybe the darker surrounding, which is free of particles, could be used as a seed for some denoising algorithm?

Manual segmentation of some samples to use as training data for a machine learning approach does not seem feasible, as it is hard to tell what's an actual particle and what is noise. That being sad, a segmentation result that looks pretty believable would be okay for now. (However, I need to do particle tracking based on such data later...)

The best I could do yet was using h-dome transform with a bit of morphological operations (it's an enlarged crop of another sample):

What would be a reasonable approach to take the long shapes and divide them into several ones? They seem to be slightly overlapping particles in reality.

Any ideas on where to go from here? Some method I don't know of?

Thanks so much!


texts:
I'm trying to classify states with a 1D-CNN-Network structure. Therefore I need to analyze the incoming vibrational data in a certain window. Two examples of such signals are shown below.

I'm continuously streaming the vibrational data. It should be possible to detect a trigger to determine the starting point of the process. This ensures that only interesting states are analyzed and classified.
My first approach for such an trigger detection was a rolling window with an average value. If this average value exceeds a certain threshold value, the trigger is set. But this window must be relatively large, because some signals stand still for a longer time (see second picture).

Each process have a different duration between 45s-85s. Therefore it would be perfect to get a trigger for the starting and for the end point of the process. The signal value for stall can be seen at the beginning and at the end of each graph. A single peak (caused from a collision) should not be detected as a trigger.

This signal is equal to the first class.

And here is a different signal for a second class.

Which other algorithms would be useful for this application?


texts: In MATLAB, i compared elapsed time to invert a Hermitian matrix using inverse(), svd(), and chol(). svd() took the longest. So is there any reason to prefer svd() to the other two methods?


texts:
I want to take a real signal on my FPGA and create a complex one.  I have FIR core implemted for a Hilbert Transform, and Xilinx's datasheet on it shows coefficients of: (-819,0,-1365,0,-4096,0,4096,0,1365,0,819) as an example.

I wanted to add more taps, so as practice I used the matlab fdesign function:

d=fdesign.hilbert(11,0.1);
hd=design(d,'firls');
coef=hd.Numerator

The coefficients seem reasonable (though I will need to multiply by 10000 to get integers).  What I don't understand is the 0.1 transition width variable (I took the 0.1 from an example online), nor if I need to tune this filter for a frequency range, or if they are generic enough.


texts: Discrete Fourier transform of even function

texts:
I have a cosine function - Hence it is even. Considering only the real parts of the DFT, on performing DFT, I am getting something like this. Could anybody tell me where I am going wrong:

DFT:

Please let me where I am going wrong with the implementation


texts: Power spectral density of surface roughness

texts: I am currently working on understanding how visually and intuitively the power spectral density is the fourier transform of the autocorrelation of a signal. I was wondering if anyone had a good explanation that would make sense to someone not directly in the field of signal processing but is now encountering it as a field I have to mesh in for a problem I'm encountering in my field.


texts: Recommended Resources / Literature Search Terms for a Solutions to a Specific Kind of Multi Harmonic Signal Structure

texts:
Hopefully this isn't considered too off-topic. I'm working in industry these days and came up with a solution to a signal processing problem we'd been facing. I'd like to get a sense as to whether said solution has been published before or if there are alternatives I should look at, but I'm having trouble searching for this particular signal scenario. I'll describe the signal structure below, and would appreciate any input on how I should be searching the literature.
So I have two signals: one that is pure noise (noise_only), and one that is very similar noise plus a target signal of interest (noise_plus_target). In each signal, the noise is actually made up of multiple somewhat-frequency-distinct signals, and the same set of said noise signals contribute to each observed signals additively but with different weights between the two. In theory, we should be able to use information from the noise_only signal to help remove the noise from the noise_plus_target signal, and I've come up with a method to achieve this, but I want to know what other solutions might have been already published for this kind of scenario.
In R, code for generating fake signals that match the characteristics of my real signals would be:
library(tidyverse)

# define a function to generate simple sinusoid given time and hz
sine = function(time,hz) sin(time*(2*pi)*hz)

#define a function to scale values to 0:1
scale01 = function(x) (x - min(x)) / diff(range(x))

#specify sample rate
sample_rate = 10 #in Hz
max_time = 30

#construct a tibble
latent_signals = tibble(
    #specify sampling times (in seconds)
      time = seq(0,max_time,1/sample_rate) #30s of data
    #construct some latent noise signals, each at a decently separated Hz
    , noise1 = sine(time,1/11)
    , noise2 = sine(time,1/3)
    , noise3 = sine(time,1)
    #specify a target signal that will be hidden in the noise
    # This could take any shape; here I've chosen a bump midway
    # through the timeseries
    , target = scale01(dnorm(time,mean=max_time/2,sd=3))
)

#show the latent signals
latent_signals %>%
    tidyr::pivot_longer(
        cols = -time
    ) %>%
    ggplot()+
    facet_grid(
        name ~ .
    )+
    geom_line(
        mapping = aes(
            x = time
            , y = value
        )
    )

#combine the latent signals into two observed signals, with different weights
# for each and the latent target only in one
latent_signals %>%
    dplyr::mutate(
        noise_only =
            noise1*runif(1,.5,1.5) +
            noise2*runif(1,.5,1.5) +
            noise3*runif(1,.5,1.5)
        , noise_plus_target =
            noise1*runif(1,.5,1.5) +
            noise2*runif(1,.5,1.5) +
            noise3*runif(1,.5,1.5) +
            target
    ) %>%
    dplyr::select(
        time
        , contains('_')
    ) ->
    observed_signals

#show the observed signals
observed_signals %>%
    tidyr::pivot_longer(
        cols = -time
    ) %>%
    ggplot()+
    facet_grid(
        name ~ .
    )+
    geom_line(
        mapping = aes(
            x = time
            , y = value
        )
    )

```


texts: Calculating penalty to SNR for non-matched filter

texts: If a matched filter is optimal for detecting a known pulse in an AWGN channel, any other filter must have a degraded SNR. Given a discrete pulse shape and an arbitrary non-matched filter, is it possible to calculate how much worse the SNR will be, relative to the matched filter case?


texts: How do I sum the discrete cosine transform cosine waves?

texts:
Sorry if this sounds like a stupid question, I'm new to DCT.

I've implemented various types of DCT (like IV and II) and they all give vastly different coefficients. From what I've read, the DCT expresses a signal as a sum of cosine waves with various oscillating frequencies. How can different formulas that give different coefficients express the same signal? I created 7 cosine waves (With one DC coefficient) on a graph with proportional periods according to N, and summed them but end up with a function that doesn't represent the signal at all. I'm probably doing this wrong, but here's what I did:

Calculated the 1D-DCT-II of this array:
 [1 2 3 4 5 6 7 8]
DCT:

[12.727922061358 -6.4423230227051 -1.3322676295502E-15 -0.67345480090394 0 -0.200902903736 -1.5765166949677E-14 -0.050702322759621]

Now I don't know if this is how to do it but I tried this (with a window of length 8):

12.72 + -6.44cos((pi/7)x)+ 0cos((pi/3.5)x)...

I can use the inverse DCT (DCT-III) and get the input values back, but I want to see the signal that represents the original input values. I tried summing the frequencies but I don't know if I'm doing it right so can someone guide me?
Btw the formula I used was

$y[$n]= $realarray[$n]*cos((M_PI/$N)*($n+0.5)*($k));


texts:
I am trying to port a MATLAB code into C++ and found the usage of rand in MATLAB.

Is MATLAB's rand function and C++ rand() function same? If not, is there any inbuilt function for MATLAB's rand in C++ or OpenCV?


texts: Suppose I have a sequence of identical symbols carrying the same information at a sample rate equal to the symbol rate, why do we always talk about averaging them to average out noise, why not just  adding them all together without  averaging them? I suppose the improvement in SNR is going to be the same assuming that noise is uncorrelated from one sample or symbol to another, is this related to the limited dynamic range in DSP?


texts:
From the Transmitter comes a block of size N+L (where N is the number of subcarriers and L is the cyclic prefix length). Then the block goes through a Channel which adds delay and AWGN, so this process can be described as a convolution of a transmitted block with a channel (the channel length H).

The question: The output of the convolution increases the size of transmitted block (the output from a channel is of lenth [N+L+H]) so how can FFT be performed (of course after removing cyclic prefix)? (because the number of subcarriers must be the size of IFFT/FFT)


texts: Inverse mapping (affine transformation) for singular transform matrix

texts:
I was reading DIP by Gonzalez et al, and came across the concept of getting output image pixel values by calculating nearest neighbour in input image by inverse mapping of output pixel coordinates $(x,y)$ to input coordinates $(v,w)$ and applying its value. Something like

$$
(v,w) = T^{-1} (x,y)
$$
where $T$ is the image transform matrix.

But what if the matrix T is singular. How will we define the nearest neighbour pixels then?
For eg. for

$$
T = \left [
\begin{array}{ccc}
1 &1 &0\\
1 &1 &0\\
0 &0 &1 \end{array}
\right],
$$ whose inverse doesn't exist.

(I am guessing that the transform matrix is changed by adding or subtracting some values)


texts: Estimation of Amplitude, Frequency and Phase of Linear Combination of Harmonic Signal Beyond the Leakage Resolution of DFT

texts:
How can I find a rough ( as accurate as possible) Amplitude of each frequency when there is spectral leakage.
Currently, I am dealing with a system that contains special leakage which seems unavoidable as I am measuring a real signal with any possible frequency. currently, I can identify the component frequencies to a decent level of accuracy even with the leakage, but I can't find the amplitude of this frequency as the energy is spread over multiple frequency bins.
Is it possible to get a rough prediction of a frequency amplitude either directly from the FFT, or afterwords using the identified frequencies?
I have been experimenting with Windowing and found this useful.
https://www.youtube.com/watch?v=VxTx9QW8Zx8&ab_channel=Adash
using the Hann window and using a frequency correction equation. it gives good approximations for frequency but no amplitude or phase.
A sample of the code currently used to test and generate the spectra.
import numpy as np

def gendata(FreqR,AmplR ,Ssize):
    PhaseR = np.pi # phase range
    FreqR = (FreqR*10) +1
    AmplR = (AmplR*100)+1
    PhaseR = (PhaseR*100)+1

    freq_OR = []
    ampl_OR= []
    phase_OR =[]

    for i in range(Ssize):
        freq_OR. append((((FreqR-1)/Ssize)*(i+1)/10)) # currently just equaly spread but can be any value between 0-9 rad/s
        ampl_OR.append(1) # unit amplitude to check for variation and error.
        phase_OR.append(np.random.randint((-PhaseR),(PhaseR))/100) # random phase

    return(freq_OR,ampl_OR,phase_OR)

def pltsignal(freq_OR,ampl_OR,phase_OR,t_OR):

    z = 0 # original signal flat sea's
    if len(freq_OR) == len(ampl_OR):
        for i in range(len(freq_OR)): # for each frequency generate a regular wave.
            wave_OR = ampl_OR[i] * np.cos(2*np.pi*freq_OR[i]*t_OR + phase_OR[i]) # generated wave signal using the frequncy, amplitude and phase of each wave.

            z = z + wave_OR # superpostion each wave ontop of each other.
    else:
        print("amplutide and frequency are diffrent lengths")

    return z  # returns superimposed irregular wave.

FreqR = 1.5 # htz or (9 rad/s ish)
AmplR = 1 # currently unit amplitude to measure the error
Size = 12 # sample size, also used to move around frequencies

E_OR = 20 # time over wich signal is measured (max values is about 120 ish)
Fs = 6 # Sample rate in Hz

t_OR = np.arange(0,E_OR,1/Fs) # starts at 0, ends at E, steps by 1/Fs

data_OR = gendata(FreqR,AmplR,Size)
z = pltsignal(data_OR[0],data_OR[1],data_OR[2],t_OR) # original signal ( simulation of mesurment )

# anaysisng the signal using NumPy fft, then scaling and finding peaks.
Ramp = np.fft.fft(z) #real aplitudes used to find phase
Rfeq = np.fft.fftfreq(z.shape[-1]) #real frequency domain to find phase

Signal Model:
$$ x \left( t \right) = \sum_{i = 1}^{M} {a}_{i} \cos \left( 2 \pi {f}_{i} t + {\phi}_{i} \right) + n \left( t \right) $$
Where $ M, {\left\{ {a}_{i} \right\}}_{i = 1}^{M}, {\left\{ {f}_{i} \right\}}_{i = 1}^{M}, {\left\{ {\phi}_{i} \right\}}_{i = 1}^{M} $ are unknown parameters and $ n \left( t \right) $ is Additive Gaussian White Noise (AWGN).
One could assume:

The SNR is very very high.
The observation time is ~ 120 [Sec].
The number of signals is 1-20.
The frequencies are up to 2 [Hz].
The gap between frequencies can be as small as 0.005 [Hz].


texts:
I'm reading a paper about texture synthesis using pyramids, and the very first section mention two methods for constructing pyramids, one is the Laplacian Pyramid (which I'm quite familiar with), the other one is "steerable pyramid". I got the idea of what the latter cope with in comparison with the former, but I do struggle to understand how they work (if I wanted to implement them from scratch for example).

Is there a sufficient simple explanation of how they work?

There's a diagram extrapolated from one of the papers I've been reading through, the intuition I have behind it is that I need to perform these operation in the fourier domain, but I might be wrong. Can anyone expound it?

Thank you


texts: What is the most precise frequency analysis spectrography algorythm?

texts:
Everyone uses Fast Fourier Transform, which is fast at the detriment of precision. The input audio has sample accuracy and the FFT has 1/64 sample accuracy.

What algorithms can output high resolution frequency analysis of audio that varies over single samples, and that can divide input audio into i.e. 1024 bands of different frequency sines, each of the 1024 tracks is 22k sample rate?


texts:
I know this question has been asked here numerous times. But I am still unclear here :(

My book says,

  A linear time-invariant system is stable if its impulse response is absolutely summable.  (G. Proakis)

So that should mean that, for a stable system, the impulse response h[n] goes to zero as n approaches infinity, right? But that filter is FIR filter. So FIR filters are always stable.

So, for unstable filters, the impulse response is not absolutely summable. In another way, the impulse response never approaches zero. Again, for IIR filter, h[h] continues to go on with n i.e. never goes to zero. So, IIR filters are supposed to be unstable.

But from Z-transform, IIR filters have poles and FIR filters don't, correct? And if the poles are inside unit circle, i.e. if their value is below one, then that filter is stable.

So, what's happening here? If IIR filters have impulse response which is not absolutely summable, how come they are stable when their poles are within unit circle? Their impulse response goes on like forever, right? Their output gets feedback. So why they can be stable depending on pole's location or value?

Thanks in advance :)


texts:
I am currently reading Candes et. al.'s 2006 paper[1] on recovery of sparse signals from incomplete frequency samples. I am having trouble figuring out what is the form of the Fourier transform restricted to $\Omega$, $\hat{f}|_{\Omega}$, which appears in Lemma 1.2.

This restriction appears to be a bijection $\mathcal{F}_{T \to \Omega}: l_2(T) \to l_2(\Omega)$, when $|T| = |\Omega|$ and $T, \Omega$ are subsets of $\mathbb{Z}_p$, $p$ prime. Here, $l_2(C)$ denotes the space of signals with $C$ as their support.

Is there a concrete example of such a $\mathcal{F}_{T \to \Omega}$, or are we only able to prove its existence? In Terence Tao's paper[2] containing the proof of this Lemma, the proof mentions that the matrix of $\mathcal{F}_{T \to \Omega}$'s coefficients has a form similar to that of a Vandermonde matrix, i.e. it will contain elements $(e^{- 2 \pi x_i \xi_j / p})$, which makes me suspect that we can obtain concrete examples of this linear transformation.


texts: Classifying and filtering out noise from zero crossings

texts:
Note: I know absolutely nothing about zero crossings nor signal processing in general. However, I am doing programming exercises that require them. I've searched around but I cannot understand most of the information I'm getting. Hopefully I get clarification here. Feel free to close if it is off topic.

Here are two exercises from a computer science book (here https://mitpress.mit.edu/sites/default/files/sicp/full-text/book/book-Z-H-24.html#%_sec_3.5.3). I'll only include the relevant details.

  Exercise 3.74: Alyssa P. Hacker is designing a system to process
  signals coming from physical sensors. One important feature she wishes
  to produce is a signal that describes the zero crossings of the input
  signal. That is, the resulting signal should be 1 whenever the input
  signal changes from negative to positive, -1 whenever the input signal
  changes from positive to negative, and 0 otherwise. (Assume that the
  sign of a 0 input is positive.) For example, a typical input signal
  with its associated zero-crossing signal would be   … 1 2 1.5 1 0.5
  -0.1 -2 -3 -2 -0.5 0.2 3 4 … … 0 0  0  0  0   -1   0  0  0   0   1  0 0 …  ...

The associated second exercise is:

  Exercise 3.75.  Unfortunately, Alyssa's zero-crossing detector in
  exercise 3.74 proves to be insufficient, because the noisy signal from
  the sensor leads to spurious zero crossings. Lem E. Tweakit, a
  hardware specialist, suggests that Alyssa smooth the signal to filter
  out the noise before extracting the zero crossings. Alyssa takes his
  advice and decides to extract the zero crossings from the signal
  constructed by averaging each value of the sense data with the
  previous value. She explains the problem to her assistant, Louis
  Reasoner, who attempts to implement the idea, altering Alyssa's
  program as follows:

(define (make-zero-crossings input-stream last-value)
   (let ((avpt (/ (+ (stream-car input-stream) last-value) 2)))
     (cons-stream (sign-change-detector avpt last-value)
        (make-zero-crossings (stream-cdr input-stream)
                               avpt))))

  This does not correctly implement Alyssa's plan. Find the bug that
  Louis has installed and fix it without changing the structure of the
  program. (Hint: You will need to increase the number of arguments to
  make-zero-crossings.)

note: sense-data is the input signal

My questions are:

What is Exercise 3.75 referring to? Is it the large amounts of 0s?
Suppose we filter out the noise from the signals in "… 0 0  0  0  0   -1   0  0  0   0   1  0 0 …", should the resulting filtered signal be "… 0 -1   0  1  0 …" or "… -1   1  …"?
In this exercise, smoothing is done by averaging each value of the input signal with the previous value, for the life of me, I don't understand what that should do and why we do it?


texts: Global partition of image sets using neighbourhood measures

texts:
I have the following general problem,

I presuppose an image is the addition of $K$ other images

$$f(x) = \sum_{k=1}^{K} f_k(x)$$

but I don't know what they are.

At each point in an image I obtain $K$ values that represent each of these images, but the order is unknown. In my particular case, these are orientation values.

$$\mathbf{g} = \{ g_j \mid j = 1 ... K\}$$

At each point, I wish to uniquely assign each value $g_j$ to each of $K$ images, based on the similarity to their neighbours. For example, for a 3x3 image with $K=2$ and $\mathbf{g}$ given by

[0,1]  [0,1]  [1,0]
[1,0]  [1,0]  [1,0]
[0,1]  [1,0]  [0,1]

the two minimal energy sets would be, for some energy function that compares differences,

0 0 0        1 1 1
0 0 0   and  1 1 1
0 0 0        1 1 1

I imagine this would require some kind of iterative optimisation process, but before reinventing the wheel I wonder if any one has come across this problem before and knows a good solution?


texts: How to downsample my data readings from 0.5hz to 0.001667hz. Do I filter and downsample?

texts: I have a device that, in the lab, takes voltage readings every 2 seconds. However in real-life application the device would only wake every 10 minutes and then take a 100 readings whilst awake. How do i change my collected data to mirror real-life application? Only then can I interpret if my device is working as it should. I see from reading various posts that downsampling alone is not appropriate.


texts: Why is White Noise so important in System Identification or Adaptive Filters

texts:
I'm looking to implement a feedback cancellation filter using Wiener Filtering, where an adaptive Wiener filter is used to cancel the feedback occurring in the path between a loudspeaker and a mic (assume PA system). The idea is essentially from this paper:

  Spriet, Ann, et al. "Adaptive feedback cancellation in hearing aids
  with linear prediction of the desired signal." IEEE Transactions on
  signal processing 53.10 (2005): 3749-3763.

According to the paper:

Signal delivered at the mic is $x[k]$.
Signal driving the loudspeaker is $u[k]$.
The transfer function of the feedback path is $F[q]$, so the signal captured at the microphone is $y[k] = x[k] + u[k]F[q]$.
The feedforward transfer function from microphone to the loudspeaker is $G[q]$.

Here is the block diagram that helps put it all together:

In the description, $k$ is the discrete time index, and $q^{-1}$ is the unit delay operator (I know these to conventionally be $n$ and $z^{-1}$).

The idea is to introduce an adaptive filter $\hat{F}[q]$ that estimates the feedback path and cancels it from $y[k]$. There's a couple of things being talked about, but essentially since $x[k]$ and $u[k]$ are correlated, they talk of adding a probe signal $r[k]$ to the loudspeaker input $u[k]$, which helps identify $\hat{F}[q]$. $r[k]$ is apparently usually a noise, which I assume helps decorrelate $x[k]$ and $u[k]$.

  My first problem is: Wouldn't the noise also be output by the
  loudspeaker? Or is it added to a copy of the signal, not affecting what is fed to the loudspeaker?

Alternately, the paper also says that many audio signals can be closely approximated as a low order AR process:

$x[n] = h[n]*w[n]$, where $w[n]$ is white noise. This condition is not satisfied for voiced speech or music, in which case a pulse train is suggested.

  So my bigger question is, why is this noise (white noise) so important
  in AR modeling, or Adaptive Filtering? It seems to defeat the purpose
  to add noise to the signal.

Any


texts: finding the audio pulse compression peak in extra-low SNR scenario?

texts: I just want to obtain some general ideas for such a scenario: if the environmental noise is extremely strong and the audio signal at the receiver is extremely weak(for instance, low range propagation), so even after pulse compression(cross correlation using linear chirp signal), the correlation peak is hidden in the noise. Are there any ways to detect the hidden peak?


texts: Why does an 8 point DFT behave differently from a 16 point?

texts:
I have a 50 Hz sine wave sampled at 1600 Samples/second.I'm computing the 16 point DFT.So my fundamental frequency is 100 Hz.

I then decimated the 1600 samples to 400 samples.I then computed the 8 point DFT using this.Now my fundamental component is 50 Hz.The second harmonic is 100 Hz.

I plotted the above two scenarios on matlab.In the first case (16 Point) the magnitude of the 100 Hz component is nearly 3.5.

In the second case (8 Point) the magnitude of the 100 Hz component is nearly zero.The magnitude of the 50 Hz component is nearly 4.

Why does this anomaly occur?

Thanks,

Jose

Here's the code for this

%----------------------------Sine Wave------------------------------%

f = 50;

t = 0:(1/1600):(1-(1/1600))

Signal = sin(2*pi*f*t);

subplot(3,1,1)

plot(t,Signal)

ylabel('Amplitude')

title('50 Hz Sine wave at 1600 samples/sec');

%-------------------16 Point DFT-----------------------%

DFT_app_16 = fft(Signal,16);

subplot(3,1,2)
stem((0:100:1500),abs(DFT_app_16))

ylabel('Magnitude')
title('16 Point DFT')

%---------------------8 Point DFT--------------------------%

Dec_Sample = [];

for i = 1:4:1597
    Dec_Sample = [Dec_Sample ((Signal(1,i) + Signal(1,i+1) + Signal(1,i+2) + Signal(1,i+3))/4)]

end

DFT_app_8 = fft(Dec_Sample,8);

subplot(3,1,3)
stem((0:50:350),abs(DFT_app_8))

xlabel('Frequency(Hz)')
ylabel('Magnitude')

title('8 Point DFT')


texts:
)Following a smoothing filter of a biological signal.
y[n] = Ay[n − 1] + Bx[n − 2]
Please find the transfer function in the z-domain and find the impulse response function of this
filter.


texts: literature related to Spectral resolution

texts:
I  am interested in spectral resolution of spectral estimation methods.

Can you point me to some good literature, videos or books from which i can read and study?

I want to understand how to resolve  two closed spaced sinusoids:
$$
x_1(t) = \sin(\omega_1 t + \phi_1) \\
x_2(t) = \sin(\omega_2 t + \phi_2) \\
$$
where $\omega_1 \approx \omega_2$.


texts:
What I don't understand is the connection between the sample rate, the baud rate, decimation and using convolution or correlation.

I'm sampling at 2048000 Hz but I don't know my baud rate, I estimated it (by looking at the time length of one of the short pulses) to be about 2000 pulses per second.

So these are my steps (roughly):

get samples
demodulate
decimate by about 1000
correlate with ones

Here is an example of what my signal looks like after demodulation:

You can see that there are 4 symbols but many many samples.  How can I turn those many many sample into just 4 symbols? I think it is related to convolution/correlation, but I'm not sure how.


texts: How to find the Fourier transfer function of a Laplacian operation?

texts:
My practice exam question is:

  A discrete Laplacian operation on an image $f(x,y)$ results in the following output image:
  $$g(x,y) = 4f(x,y) - f(x+1,y)-f(x-1,y)-f(x,y+1)-f(x,y-1)$$

  Find the 2 dimensional discrete Fourier transfer function of the Laplacian operation and show that the laplacian operation is a high pass filter.

I have the answer sheet, but I really don't understand it so am asking for a bit more of an explanation on the steps to solve it.


texts: why is my FFT magnitude of real time recording so volatile?

texts:
My goal is to do spectrum analysis of a gear motor. I used Python with the Pyaudio package and made a crude spectrum analyzer that displays the FFT as the sound is being recorded. I recorded 10 seconds of sound just for test. Sample rate is 44100, each frame is 1024 samples. For each frame I plotted the FFT, and made a dynamic plot thru the 10 seconds.

I also downloaded a FFT analyzer to my smartphone and compared the results as I turned on a motor. The smartphone app is able to display relatively stable spectrum. In my plot, the magnitudes are very high upon turning on the motor, and only "settled" after a couple of seconds.

I wonder why that is. Is it due to my sound card or laptop mic? I tested using another laptop and got similar result. Is it because I recorded in Mono mode? The spectrum changes too much as I move the motor around, whereas on the smartphone the spectrum is more stable.

Or is it because I didn't use a Window function? Here is my code:

import numpy as np
import os
import time
import pyaudio
import matplotlib.pyplot as plt
from scipy.fftpack import fft

CHUNK = 1024
FORMAT = pyaudio.paInt16
CHANNELS = 1
RATE = 44100
RECORD_SECONDS = 10
WAVE_OUTPUT_FILENAME = "test.wav"

p = pyaudio.PyAudio()

stream = p.open(format=FORMAT,
                channels=CHANNELS,
                rate=RATE,
                input=True,
                frames_per_buffer=CHUNK)

input('Press key to start recording:\n')
print('1 sec delay started') #delay so keystroke doesn't get recorded
time.sleep(1)
print("* recording")

frames = []
n = 1024
k=np.arange(n)
T = n/RATE
frq = k/T
frq = frq[range(int(n/2))] # one side frequency range
for i in range(0, int(RATE / CHUNK * RECORD_SECONDS)):
    data = stream.read(CHUNK)
    frames.append(data)
    decoded = np.fromstring(data, dtype=np.int16) #grab the data in stream
    fft_decode=fft(decoded)/(len(decoded)/2) #normalized FFT
    mags=np.absolute(fft(decoded)) #
    plt.ylim(top=55000)
    plt.xlabel('Freq (Hz)')
    plt.ylabel('|Y(freq)|')
    plt.plot(frq, mags[range(int(n/2))],'b')
    plt.pause(.01)
    plt.gcf().clear()

print("* done recording")

plt.close()

How can I improve my analyzer?
 Is there anything wrong with my code or the way I approached this?

Edit: Nevermind. The reason turns out to be the boost function of my microphone, and disabling it solved my problem.


texts: How does Adobe After Effects generate its "audio spectrum" effect?

texts:
I'm trying to replicate the "audio spectrum" effect from Adobe After Effects. An example can be seen in this video:

Obviously, it has to be some variant of a fourier transform, but I've tried taking the fft of an audio sample in realtime and my implementation looks more like a big block of bars:

I imagine some kind of intermediate function is being applied to the fft, but I can't really figure out what it might be.


texts:
So I've been looking at this butterfly diagram to try to understand it better:

And I am trying to get a good understanding of the twiddle factors.
The definition is given as:
FFT Twiddle Factor: ${e^{i2{\pi}k/N}}$
and
IFFT Twiddle Factor: ${e^{-i2{\pi}k/N}}$
So k is the index number of the iteration thus $k=0,1...N$ but its $N$ that I am unsure of.
From the image is the first stage N = 8 (since there are 8 butterflies) or is N = 2 since each butterfly only spans two elements? Or is N always 16 every pass?
My guess currently is:
Stage 0 : N = 8
Stage 1 : N = 4
Stage 2 : N = 2
Stage 3 : N = 1

Am i correct in this assumption? If not i hope some one clarifying my misunderstanding here.
Thanks


texts: Simulating noise in computed tomography reconstruction

texts:
In the research of the computed tomography (CT) reconstruction, one needs to simulate the noise during CT projection and capturing. Then, my questions are:

How about the proper noise type?
How to add the noise, in the sinogram domain or the frequency domain?

Any reference materials are also preferred.


texts:
Problem: I am looking at an adaptive filtering application where the eigenvaluespread of the autocorrelation matrix $R$ is important for the convergence of the algorithm. For a single channel system the autocorrelation matrix $R$ for iterationstep $n$ can be calculated by $R=E\{ x(n)  x^H(n)\}$ where $x(n)$ is the input signal of the adaptive filter at iterationstep $n$ consisting of a number of samples $N$ recorded over a timespan. The calculation of the eigenvalues is straight forward.
Question: What is the "multichannel equivalent" for $R$ in the case of e.g. an adaptive filtering multichannel application? Do I need to calculate some sort of autocorrelation tensor?


texts: What does Fast and Slow Time Weighting signify for real-time audio signal processing?

texts:
I have a streaming audio signal that I am recording in real time. To make the explanation easy, say the signal that I receive is calibrated and has a sampling frequency, $F_s = 44100\textrm{ Hz}$.

I am then finding the Sound Pressure Level of this signal. For this purpose, I get a frame of audio data of t seconds.

Then, RMS Sound Pressure Level (dB) is given by:
$$
\rm Amplitude_{RMS} = \sqrt{mean\left(frame ^2\right)}
$$

This gives the instantaneous RMS of a $t$-second frame of audio. Then, to find a $\textrm{ dB}$ level,
$$
L_p = 20\log_{10}\left(\rm Amplitude_{RMS}\right)
$$

Things get tricky here with the definitions of 'Fast' or 'Slow' sound pressure level. As this post points out, for time-weighted exponential average, I created a low-pass filter with a real pole at $1/\tau$ (where $\tau= 125\textrm{ ms}$ for 'fast' or $1\textrm{ s}$ for 'slow')

τ = 0.125; %fast
[b a] = bilinear(1, [1 1/τ], Fs);
frame_filtered = filter(b, a, frame.^2);
L_fast = 20*log10(sqrt(mean(frame_filtered)));

The IEC 61672 specifies the time weightings to be Fast and Slow for determining the speed at which the instrument responds to changing noise levels.

How does this work for digital audio?
Does it mean that the analysis needs to be performed on a fixed frame of data?
For this particular example, this frame would be $t\cdot F_s$ samples?
Slow: 1*44100 = 441000 samples
Fast: 0.125*44100 ~= 5513 samples
If so, is it true that I would need to know in advance whether I need fast or slow levels?


texts:
I need to compare images using Gabor filter. I cant seem to understand how to find a quantifiable value to compare through Gabor.
Can someone give me steps to do the same. I got some code, but it doesnt compare two images.
Any help is appreciated.


texts: inter symbol interference vs digital modulation

texts:
I am a little confused  after my lecture today. My professor started off with explaining what inter symbol interference (ISI) is and he ended the lecture talking about the different types of digital modulations out there.

I am now confused with how these two actually come together, or is it two different topics that have no connection?

Is ISI about sending digital signals in a wire? Does ISI come into play when doing lets say QPSK? I thought the symbol signal sent in QPSK is just a sin signal. or it uses this ISI?


texts: This article mentions that the DFT/FFT can be used as a channeliser in a similar way to a filterbank. I get how the DFT filters into spectral bins, but its output is in the frequency domain, i.e. Fourier coefficients, nolonger in the time domain, so either I'm missing something or they've neglected to mention that DFT filterbank output would need to be transformed back to the time domain to function like an actual filterbank. Can someone shed light on this?


texts: object tracking in raw h.264 compressed video

texts:
I am working on a project that requires me to detect and track a human in a live video from a webcam connected to a Beagleboard xm. I have completed this task using Opencv in pixel domain. The results on the board are very accurate but extremely slow. Many people have suggested me to leave pixel domain and do the same task in an h.264/MPEG-4 compressed video as it would extremely reduce the computational overhead. I have read many research papers but failed to discover any software platform or a library that I can use to analyze and process h.264 compressed videos. I will be thankful if someone can suggest me some library for h.264 compressed video analysis and guide me further.
I know about P and I frames but don't know how to extract them and use which programming software or language to process them.

Thanks and Regards.


texts: Can disparity range from negative to positive values for parallel cameras?

texts:
In my stereo-rectified camera pair (not physical, but two positions of the same camera), when I project some 3D points in both images to find a rough disparity range, I often get a range that goes from negative to positive (e.g. -11, 24). I have no problem with ranges like (-22, -10) or (10, 15).

Is that kind of range even possible? When I do stereo matching how should I treat negative values? And what about the depth calculated from that disparity that can have negative values (and possibly areas with lowest cost at 0 disparity)?


texts:
I am studying the Kalman filter and its basic implementation, and it was asked to use the filter to estimate a signal observed in noise
$$y(n) = x(n) + v(n)$$
where $v(n) \sim \mathcal{N}(0, \sigma^2)$ and $x(n)$ is modeled as an ARMA(2,2) process
$$x(n) = b_0 u(n) + b_1 u(n-1) - a_1 x(n-1) - a_2 x(n-2).$$

The exercise is conceptually intersting, but I was wondering if there is an direct application for this case. I mean, which processes, in real life, can be modeled as an ARMA(2,2) and in which context the Kalman filter is used, not as a predictor, for a signal like $s(n)$?

Thanks!


texts: Tricks for plotting the magnitude of a DTFT?

texts:
To preface, this isn't a homework question but rather a self-study question to help me to understand the basics of finding the DTFT and magnitude of the DTFT based on a discrete time signal sampled from a continuous time signal.

I was wondering if anyone has any methodology to approaching questions where you're asked to plot the magnitude of the DTFT of a discrete-time signal, when
$x(t)$ is sampled according to $x[n] = x(nT_s)$`, with $F_s$ equal to a given value.

For example if $$1/T_{s}=8W$$ and :
$$x(t) = T_{s}\frac{1}{4W}\frac{\sin(2\pi Wt)}{\pi t}\frac{\sin(2\pi 2Wt)}{\pi t}$$

What are the first questions that you would ask yourself that will help you to plot the magnitude of the DTFT.  I am having trouble translating the DTFT to the plots themselves.


texts: How to identify causality, stability and ROC from the pole-zero plot?

texts:
To preface, this is not a homework related question but purely for self-study purposes.

If I am given the following Pole-Zero Plot: (Source: Berkeley Exam1)

How would I go about trying to determine whether the system is causal or stable?
Similarly, why should the region of convergence be $0.5 < |z| < 2 $ for stability?

I know that for stability, the ROC must include the unit circle.  And for causality, $h[n]=0$ for $ n<0$.  Other than that, I am having trouble applying those rules to a pole-zero plot.


texts: Effect of DC component on the whole signal - comparison between normalised and non normalised

texts:
I have a fourier analysis signal as in the picture attached, where red represents the FFT of movement of the hand of a stroke subject and the blue one is the movement of a healthy subject.

I am doing some analysis called Spectral Arc Length, where I will calculate the spectral arc length to compute the smootheness of movement. (check out this paper: On the analysis of movement smoothness by Sivakumar Balasubramanian), where in the metric, the longer the spectral arc length, the less smooth the movement is.

The first image down here is the original data, where we can see that the DC component of the stroke subject is higher than the healthy subject, and the amplitude of the frequency signal is higher, causing the arc length of the signal to be larger, signifies less smooth movement.

MATLAB CODE:

y=variable;
Ts=1/Fs;
L= size(y,2); % Length of signal
NFFT = 2^(ceil(log2(L))+4); % Next power of 2 from length of y
Y = fft( y, NFFT )/L;

f = Fs/2*linspace(0,1,NFFT/2+1);

however, I have a problem with my signal where the metric suggest than we should normalised the signal to the DC component of the signal (Y=Y/max(Y)). By doing this, my signal turn out to be like the second picture, and when I calculate the spectral arc length metric, turns out that stroke 'apparently' have smoother movement than the healthy (due to shorter spectral arc length), which I am pretty sure should't be the case.

My question is, does the normalization to the DC component makes sense? Does the DC component has any effect on the rest of the signal, where larger DC will cause larger amplitude of the signal?
another option for me is to calculate the arc length after the DC signal (starting from the black marker is put on the signal here)
I would also try to do some wavelet analysis if Fourier analysis doesn't work for my calculation...

Non-normalised data

Normalised data Y= Y/max(Y)


texts: Understanding the magnitude of frequency response filterbank based on elliptic filters

texts:
I've implemented a 10-channels filterbank with octave-scaled second-order elliptic filters using the Python's library scipy.signal.

Here is the magnitude of the frequency response:

Can someone please explain to me why the response changes with increasing frequency? I was expecting each filter magnitude response similar to the high-frequency ones (more or less from 800 Hz upwards). Particularly, I was expecting a different response from the first 4 filters.

Here is the code:

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal

fs = 44100

# Elliptic filters
subbands = {'Band01': [50],
            'Band02': [50,100],
            'Band03': [100,200],
            'Band04': [200,400],
            'Band05': [400,800],
            'Band06': [800,1600],
            'Band07': [1600,3200],
            'Band08': [3200,6400],
            'Band09': [6400,12800],
            'Band10': [12800]}

band_ticks = [10,50,100,200,400,800,1600,3200,6400,12800,22050]
labels_ticks =["0","50","100","200","400","800","1.6k","3.2k","6.4k","12.8k","22.05k"]
filters = {}
f_resp = {}
for key in subbands:
    wc = [wc/(fs/2) for wc in subbands[key]]
    if key == 'Band01':
        filters[key] = signal.ellip(N=2,rp=3,rs=60,Wn=wc,btype='lowpass',analog=False)
    elif key == 'Band10':
        filters[key] = signal.ellip(N=2,rp=3,rs=60,Wn=wc,btype='highpass',analog=False)
    else:
        filters[key] = signal.ellip(N=2,rp=3,rs=60,Wn=wc,btype='bandpass',analog=False)
    f_resp[key] = signal.freqz(filters[key][0],filters[key][1])
    mag = 20*np.log10(abs(f_resp[key][1]))
    freqs = f_resp[key][0]*fs/(2*np.pi)
    plt.semilogx(freqs,mag,'b')
    plt.ylim([-35,4])
    plt.xlim([10,fs/2])
    plt.grid(axis='both', linestyle='-', color='grey')
    plt.xticks(ticks=band_ticks,labels=labels_ticks)
    plt.xlabel('Frequency (Hz)')
    plt.ylabel('Magnitude (dB)')


texts:
Many applications that produce sound, such as software synthesizers, are able to apply a filter that varies with time, such as applying a low pass filter that varies with an LFO.

I currently have a means of calculating coefficients for an IIR filter given poles and zeros, but do not know how to achieve this varying without redoing the calculations from scratch each time it needs to change. Is there a good approach to this that does not rely on using only a first order filter?


texts:
I'm trying to solve the following nonlinear cross-correlation problem for the time-domain signal $f(t)$:

$S(\omega) = \overline{\mathcal{F}\left[f(t)\right]} \mathcal{F}\left[f^n(t)\right]$

with $n>1$, the bar denoting the complex conjugate, and $S(\omega)$ is a given (complex-valued) spectrum. $f(t)$ is the unknown in this implicit equation and I want to solve for it.

I know the problem is undetermined for $n=1$ because the phase cancels out. However, for $n>1$, the phase is nontrivial and the solution for $f(t)$ is somewhat constrained. With a good enough initial guess I can use an iterative algorithm, which converges to the correct solution. However, I'm wondering if there's a smarter way to "decorrelate" $f(t)$ by exploiting the relationship between the two parts of the cross-correlation. Certainly, there is no unique solution (e.g. it's invariant against delay or phase shift), but maybe there's an analytical way to retrieve one of the solutions or a partial solution.

Thanks for any ideas!


texts: How to Model the Dispersion Which Is Introduced by an Electro Absorption Modulator (EML)?

texts:
I need to model the dispersion which is introduced by an EML in Matlab. The parameters are:

Lambda: 1310nm

PAM4 power:
Pmax = 1.09 mW
Pmin = 0.13 mW

Bandwidth = 56Gbps


texts: Time reversal of input signals in convolution

texts:
Let $f(t)$ and $g(t)$ be signals.

I can't figure out what happens to the convolution of $f(t)$ and $g(t)$ if

both $f(t)$ and $g(t)$ are reversed
one of $f(t)$ and $g(t)$ is reversed

I tried the change of variable in integral, but couldn't find out the answer.


texts: PortAudio: how to generate multiple sounds

texts:
I have simple wave generator in PortAudio. I'm creating single sine wave, and then I use a lookup sine to play it. Now I've stucked here. I have no idea how I can add another sine sound. Is it possible to somehow add new channel with another sound?

Here is my code:

soundengine.cpp

#include "soundengine.h"
#include <QDebug>

void SoundEngine::run()
{
    PaStreamParameters outputParameters;
    int i;
    double t;
    turnOFF = false;
    isPressed = false;

    static unsigned long n=0;
    for( i=0; i<TABLE_SIZE; i++, n++ )
    {
        t = (double)i/(double)SAMPLE_RATE;
        data.sine[i] = 0.80f * (float) sin( ((double)i/(double)TABLE_SIZE) * M_PI * 2. );
    }
    data.phase = 0;
    data.phase_increment = CalcPhaseIncrement(0);

    err = Pa_Initialize();
    if(err != paNoError) qDebug()<<"Błąd przy inicjalizacji strumienia:"<<Pa_GetErrorText(err);

    outputParameters.device = Pa_GetDefaultOutputDevice(); /* default output device */
    if (outputParameters.device == paNoDevice) qDebug()<<"Błąd: Brak domyślnego urządzenia wyjścia!";

    outputParameters.channelCount = 2;       /* stereo output */
    outputParameters.sampleFormat = paFloat32; /* 32 bit floating point output */
    outputParameters.suggestedLatency = Pa_GetDeviceInfo( outputParameters.device )->defaultLowOutputLatency;
    outputParameters.hostApiSpecificStreamInfo = NULL;

    err = Pa_OpenStream(
              &stream,
              NULL, /* no input */
              &outputParameters,
              SAMPLE_RATE,
              FRAMES_PER_BUFFER,
              paClipOff|paDitherOff,      /*paNoFlag we won't output out of range samples so don't bother clipping them */
              patestCallback,
              &data );
    if(err != paNoError) qDebug()<<"Błąd przy otwieraniu strumienia:"<<Pa_GetErrorText(err);
    err = Pa_StartStream( stream );
    if(err != paNoError) qDebug()<<"Błąd przy starcie strumienia:"<<Pa_GetErrorText(err);

    while (turnOFF == false) {
        Pa_Sleep(500);
    }

    err = Pa_StopStream( stream );
    if(err != paNoError) qDebug()<<"Błąd przy zatrzymywaniu strumienia:"<<Pa_GetErrorText(err);
    err = Pa_CloseStream( stream );
    if(err != paNoError) qDebug()<<"Błąd przy zamykaniu strumienia:"<<Pa_GetErrorText(err);
    Pa_Terminate();
}

float SoundEngine::LookupSine( paTestData *data, float phase )
{
    float fIndex = phase*TABLE_SIZE;
    int   index = (int) fIndex;
    float fract = fIndex - index;
    float lo = data->sine[index];
    float hi = data->sine[index+1];
    float val = lo + fract*(hi-lo);
    return val;
}

int SoundEngine::patestCallback(const void *inputBuffer, void *outputBuffer, unsigned long framesPerBuffer, const PaStreamCallbackTimeInfo *timeInfo, PaStreamCallbackFlags statusFlags, void *userData)
{
    paTestData *callData = (paTestData*)userData;
    float *out = (float*)outputBuffer;
    float sample;
    unsigned long i;

    (void) timeInfo; /* Prevent unused variable warnings. */
    (void) statusFlags;
    (void) inputBuffer;

    for( i=0; i<framesPerBuffer; i++ )
    {
            sample = LookupSine(callData, callData->phase);
            *out++ = sample;  /* left */
            *out++ = sample;  /* right */

        callData->phase += callData->phase_increment;
        if( callData->phase >= 1.0f ) callData->phase -= 1.0f;
    }

    return paContinue;
}

void SoundEngine::generateSine()
{
    if(isPressed == false)
    {
        data.phase_increment = CalcPhaseIncrement(220.0);
        isPressed = true;
    }

}

void SoundEngine::removeSine()
{
            data.phase_increment = CalcPhaseIncrement(0.0);
            isPressed = false;
}

To play sound I use generateSine() method. Maybe is there better way to play sound at given frequency, using one method(like generateSine() in my example). I would like to play two sounds by using generateSine() method twice.


texts: Solve for minimum and maximum value of dct coefficient


texts: Minimum number of bits required for each character

texts:
Assuming that a computer is interfaced with a device that has 210 characters and
that each character is transmitted using binary words what is the minimum
number of bits needed to represent each character? To find the minimum is there any specific formula ?


texts: How much can the FFT noise affect the timbre analysis of a harmonic sound?

texts: The timbre of an harmonic sound (e.g. the steady part of a violin note) is defined by the amplitude of each harmonic (including the fundamental). In order to compare the timbre of two instruments (e.g. two violins playing the G string), I could compare those amplitudes. If the instrument are similar, I suppose that there will be a little difference in the harmonic distribution. Now, is the FFT a reliable tool for this kind of analysis? Should I take into account the noise that it introduces, or is it negligible?


texts: Multivariate Histogram Matching for Image Color Matching

texts:
Suppose X is random variable that is spatial position in finite range, which its Cumulative Distribution Function (CDF) could be calculated easily, then if I have two random variable X and Y spatial position, can I make multivariate CDF by assuming those pair as one dimensional vector in the pair possible finite range?
TL;DR
I am working on experiment to match 2 images color hues, inspired from histogram matching by using CDF of each image, and create lookup table (LUT) from it. The naive way in each image, one channel used as the random variable (X), when x is the color range [0..255], so CDF will satisfy f(x) = P(X <= x), and to match both function is to map f_src(x) to f_dst(x).
What I am trying to do is to transform source image distribution in the 2d hue space to match target image distribution, by taking Cb and Cr component of each image, treat it as position in CbCr space (p,q), getting its CDF, and do histogram matching. I attempted to treat (p,q) as flattened points (p+q*256), that is element of 'CrCb range' [0..65536], build CDF, and find the LUT.
source = np.asarray(Image.open('./misa_2.png').resize((640,480)).convert('YCbCr'))
target = np.asarray(Image.open('./misa_3.png').resize((640,480)).convert('YCbCr'))

# source CbCr as flattened pair of points
flatpt_src = source[:,:,1:].reshape(-1,source[:,:,1:].shape[-1]).swapaxes(1,0).copy()
flatpt_trg = np.concatenate((flatpt_trg, np.array([[0,0],[255,255]]).T), axis=1)
flatpt_src.sort(axis=1)

# Cb -> x, Cr -> y
x_src_flat = (flatpt_src[1,:] + flatpt_src[0,:]*256) # how [1,255] will be like if flattened
y_src_cdf = np.arange(1, x_src_flat.shape[-1] + 1) / x_src_flat.shape[-1] # 1/len ... len/len

# target CbCr as flattened pair of points
flatpt_trg = target[:,:,1:].reshape(-1,target[:,:,1:].shape[-1]).swapaxes(1,0).copy()
flatpt_trg.sort(axis=1)

# Cb -> x, Cr -> y
x_trg_flat = (flatpt_trg[1,:] + flatpt_trg[0,:]*256) # how [1,255] will be like if flattened
y_trg_cdf = np.arange(1, x_trg_flat.shape[-1] + 1) / x_trg_flat.shape[-1] # 1/len ... len/len

# Source CbCr CDF, Target CbCr CDF
plt.plot(x_src_flat,y_src_cdf) # blue
plt.plot(x_trg_flat,y_trg_cdf) # orange

lookup table from projecting cdf source to cdf target
# check if lut plot shape is same as 'cdf_src(x) = cdf_trg(x)' mapping
plt.plot(flatpt_lut_index, lut2D) # blue
plt.plot(ka, el) # orange

when LUT is applied to source, this is the result
result = np.dstack((source.copy()[:,:,0], CbCr_res))

axshows = (source[:,:,1], source[:,:,2], result[:,:,1], result[:,:,2])
axtitles = ('source Cb', 'source Cr', 'target Cb', 'target Cr')

f, axs = plt.subplots(2,2, figsize=(12.8,7.2), dpi=100)
f.tight_layout()
for ax, ashw, attl in zip(axs.flat, axshows, axtitles):
    ax.imshow(ashw, cmap='gray')
    ax.set_title(attl)

the result is clearly what I didn't expecting, and at this point I am confused whether which is wrong, is it my 'multivariate' workaround or just the whole implementation.
(edit : 01/28 : update CDF plot and CbCr result-source comparison)


texts: Frequency domain samples to time domain without knowing the sampling rate?

texts:
I have some frequency response data from 802.11a OFDM communication in channel 8 of the 5 GHz band, and I would like to inverse transform this to produce the corresponding time domain response. That is, I have $N$ equally spaced samples $H(f_1), H(f_2), ..., H(f_N)$ (except the middle sample is missing, see below) of the radio channel transfer function, and I wish to obtain the time domain impulse response $h(t)$, or in practice a discrete approximation $\hat{h}[n] \approx h[n] = h(nT_s)$.

Now, if I inverse FFT the samples I get a sequence $g[n]$, but I do not know what time values the elements of the sequence represent. The ways I know of to obtain the time resolution are all based around knowing the sampling frequency $f_s = 1 / T_s$ that was used to produce the original FFT, but that is unknown here. I do however know the frequencies $f_1, f_2, ..., f_N$, and intuition tells me that it should be possible to infer the time resolution from this knowledge. I am not very familiar with discrete transforms however, so I haven't been able to figure it out myself.

To be specific, the $N = 52$ frequencies in question are centered around $5040 \; MHz$ and spaced $312.5 \; kHz$ apart. The center frequency is missing (unused), so there are 26 below center and 26 above, making the full range $f_1 = 5031.875 \; MHz$ to $f_{52} = 5048.125 \; MHz$.


texts:
I have been doing some research, because I used this algorithm that's very spread in the web:

y += a * (x - y)

or

y = (a * y) + (x - (a * x)

Where x is input, and y is the filter output, and a is actually alpha, which there's some math behind it to define it's value.

This algortihm is seen on wikipedia, this article and many more.

But none actually shows why this algorithm is a low pass filter and from which formula it actually came from(wikipedia has no source for it). I'm doing a small academic paper, so I need to know exactly from where it came.


texts: Relative merits of fixed point vs floating point computation?

texts:
I have a digital signal processing system that operates on a fast x86 machine using double precision floating point numbers.
It occurred to me that I am not really using the huge dynamic range of the floating point representation -- all quantities fit easily in the range ±32768.

My question: is it possible that switching to fixed point computations would provide a benefit in numerical accuracy (high priority) or computation time (low priority)?

Of course, the answer depends on how many bits are available for a fixed-point calculation.  How many bits of precision do typical fixed-point systems utilize?  Is it possible to efficiently do fixed-point calculations, with, say, 64 bits (16 bits integer part, 48 bits fractional part) on x86-64?

I had always thought that fixed-point calculations were used only in situations where CPU power is limited--does it make sense to use fixed-point calculations when CPU power is no concern?


texts: Extrapolation after taking inverse Fourier transform

texts:
Suppose I have 2000 point(time series) signal and I take 512 point FFT, after this some operation is done in frequency domain. Now if I take IFFT to get 512 points in time domain, Is there a way to map the output 512 IFFT points to the original time domain 2000 points.

Inverse  of this problem would be 400 point(time series) data and 512 point fft, and finally 512 IFFT points, here we can use Overlap and Add method to get 400 point output, is there any similar method for above problem.


texts: Good Continuing Education Course in the Basics of Frequency Analysis

texts:
All,

Are there any good continuing education courses of length 2-3 days to give an engineer a good background in Frequency Analysis.  B&K used to teach one, but I don't think that they teach it anymore.  I am looking for something that teaches sampling, aliasing, continuous vs. discrete signals, DFT, PSD, windowing, spectral leakage, filtering, etc.

Thank you.


texts: Normalize harmonics to the same fundamental frequency in a time series using MATLAB

texts:
I'm a bit new to DSP (so apologies if I get the terminology incorrect) and was wondering if the problem I have is standard. I have measured a time series of a fan and noticed that the my spectrogram contains harmonics. I haven't got a picture to hand but this picture I found demonstrates exactly what I'm seeing:

(First 2.0 seconds)

What I would like to do is normalize my data so that all the harmonics line up on the same respective frequencies: I want to enforce the same fundamental frequency on all my data.
This would tell me that the shifts in frequency are harmonically related (and therefore, can be linearly scaled) or that they're not harmonically related.

Is there such a filter / technique / methodology which does such a thing?

If there isn't then the solution I propose is to perform a peak detection on each slice of the spectrogram, and then calculate the linear coefficients from the peak positions.

Hopefully I've explained my problem succinctly.


texts:
I have looked at various sources, but still have not been able to answer the question: what is compression ratio? I have come to the understanding that this is the definitive answer:

compression ratio = (block size * bits per pixel) / (dct coefficients * bits per pixel)

--or--

compression ratio = amount of bits in original / amount of dct coefficients

Please confirm.


texts: How to Zero Pad in Order to Perform Filtering in the Fourier (Frequency) Domain?

texts:
Consider an $M\times N$ image $f$ and an $G \times K$ filter $h$. Given that convolution in the spatial domain corresponds to multiplication in the Fourier domain, then we can perform a convolution of $f$ with $h$ (that is, we can filter image $f$ using filter $h$) by multiplying the Fourier transforms of $f$ and $h$, $\hat{f}* \hat{h}$, and then perform an inverse discrete Fourier transform of the result, $\hat{f}* \hat{h}$, where $\hat{f}$ and $\hat{h}$ are respectively the Fourier transform of image $f$ and kernel $h$.

However, to do that, we need to make sure that $\hat{f}$ and $\hat{h}$ have the same dimension. So, we use the zero padding technique. We pad both the image $f$ and $h$, before finding their Fourier transforms. More precisely, we pad each of them with zeros, such that their new size is $(M + G - 1) \times (N + K - 1)$. This padding is apparently required (for both $f$ and $h$), but I am trying to understand why.

Why exactly is this the case? I think it has something to do with the assumption that $f$ is periodic. So, the conclusion would be that we can't simply zero pad $h$ (without also zero padding $f$), but I don't get why.


texts: It is possible to obtain non-zero padding fft from zero padding fft?

texts:
This is just a question out of my curiosity, normally we do not do this
I have a x(n)=[1 2 3]

y(n) = fft[1 2 3 0] =  6 + 0i  -2 - 2i   2 + 0i  -2 + 2i %zero padding
y'(n) = fft([1 2 3]) = 6 + 0i  -1.5 + 0.87i  -1.5 - 0.87i

I used Matlab to generate these results.
My question is there way to obtain y'(n) from y(n)?


texts: Calculate the intersample peak of a periodic sequence?

texts:
Similar to Calculating the PDF of a waveform from its samples, but for periodic sampled signals, and I just want the peak.  Is there a way to calculate the highest inter-sample peak of the reconstructed signal?

For instance, the signal 0, 0, 1, 1, 0, 0 is symmetrical, so you know the intersample peak is between samples 2 and 3, at 2.5, and can be calculated as ... + sinc(-0.5) + sinc(0.5) + sinc(5.5) + sinc(6.5) + ... ≈ 1.24

But is there a way to do this for the general case, when you don't know where the peak occurs?

This page says the MLS has the lowest possible crest factor, but I think for analog use, this isn't actually true, since MLS has strong intersample peaks, >6 dB higher than the samples, and a sweep-like signal like the Newman method has better crest factor.


texts:
Is there a way to design a lowpass FIR filter with a more relaxed transition band than what MATLAB FilterDesigner tool generates?
What I intend to do is to reduce number of taps needed to implement the filter. The response that MATLAB generates is good but it needs too much coefficients and the response that I need can be relaxed in the almost second half of the transition band.
The spec. of the desired filter is as below:
Fpass = 5 M;
Rejection at 700k offset from band edge: 15 dB
Rejection at 1.5M offset from band edge: 30 dB
Rejection at 5.5M offset from band edge: 70 dB
Here are the parameters I have used (Units are MHz and dB).
Fs = 187.5;
Fpass = 5;
Fstop = 10.5;
Apass = 0.35;
Astop = 70;
Other parameters are method = equiripple and density factor = 20;
Though the resulted filter has 70dB rejection at 10.5 M, I doesn't have the rejections needed at 700k and 1.5M offsets. To get that 15 dB rejection I have to either increase rejection at stopband or decrease Fstop which will result in more coefficients and is overkill for rest of the filter.
Also tried firpm() and using the parameters below I got the same response as the one by FiterDesigner which didn't help.
>> [n,fo,ao,w] = firpmord([5e6 10.5e6],[1 0],[0.01 0.0002],187.5e6);
>> b = firpm(n,fo,ao,w);


texts: How to remap center frequencies of a mel filter bank?

texts:
I'm trying to implement vocal tract length peturbation in python (Vocal Tract Length Perturbation (VTLP) improves speech recognition. Navdeep Jaitly, E. S. Hinton, Published 2013, Computer Science) by shifting the central frequencies of a mel filter bank (as described in the paper) but I can't wrap my head around how to even apply a simple change to the center frequencies. I'm using librosa to get the mel filter bank as a num_mels * (NFTT / 2 + 1) matrix but can't think about how to transform this matrix in such a way that shifts the center of the mel filters. Applying some scalar transform elementwise to the matix transforms the amplitudes (which makes sense to me) but I can't see how you'd transform the frequencies?
Any help on this would be greatly appreciated.