texts:
I've done research and could not find any resources on this topic.

I have message signal composed of more sinusoids. Then I preform amplitude modulation, that is multiply that signal with higher frequency carrier signal, and obtain modulated signal. Then to demodulate that signal, I multiply the signal with same carrier signal, and obtain demodulated version of signal. Then I should pass the signal through low pass filter, to obtain signal identical to message signal.

If I do this by matlab using fdesing.lowpass it is fine I get good output, but how to do this low pass filter by hand on paper, if I just ignore frequencies higher than the maximum frequency in message signal I don't get same amplitudes.

I multiply signals by using inverse Euler formula for $\cos$.

Any link or help would be appreciated.

If I didn't make myself clear, let me now and I will try to write an example (which I omit for now because I don't see here any latex support for writing mathematical equations).


texts:
I think i am not understanding the concept of convolution well.

Lets say we are given a system impulse response in the S-domain, and we have implemented a controller $G_c(s)$ that will adjust the system response to some kind of input. Now when we are implementing this controller say using an operational amp. circuit all we need to do is to get the differential equation of $G_c(s)$ and start implementing it with the correct component values. Is this correct?

And if so, lets say we are trying to implement this controller on a computer, where our inputs, outputs are no longer continuous signals. And our controller in the Z-domain is represented by $G_c(z)$. Now the next step to implement this controller in a computer is to get the difference equation and implement it in a software. Is this correct?

And if so, what is the difference between implementing the controller difference equation on a software and implementing the convolution algorithm on the software?

For example if we are given g[n] as an array of samples representing the the controller unit impulse response, and our input is stored inside x[n] a simple for loop might do this algorithm for us. As far as i know

for i in len(x):      #loop through the input elements
 for j in len(g):     #loop through the controller unit impulse response elements
  y[i+j] += x[i]*g[j]  #Store the conv. result in the output array y

Which is just telling the software to scale every point in the impulse response by x[n] and then shift the response by n, then add them all up.

I am probably missing something here !


texts: Impulse response of the inverse system to the backward difference

texts:
I'm currently stuck on a problem. I have a DT LTI system $\mathcal{T}$ that depicts the first backward difference:

$$y[n] = \mathcal{T} \{ {x[n] \} } = x[n] - x[n-1]$$

Its impulse response is therefore: $h[n] = \delta[n] - \delta[n-1]$, where $\delta[n]$ is the DT unit impulse.

I want to compute the impulse response $h'$ of the system $\mathcal{T}'$ that, when cascaded with $\mathcal{T}$, recovers my input:

$$\mathcal{T}'\{\mathcal{T}\{x[n]\}\} = x[n]$$

I know the answer is that $h'[n]$ is the DT unit step $u[n]$, but I can't quite figure out why.

Here is my idea so far:

I know that the impulse response of the cascaded system is: $$h_{casc}[n] = \mathcal{T}'\{\mathcal{T}\{\delta[n]\}\} = \delta[n]$$
I know that this impulse response is the convolution of the impulse responses of $\mathcal{T}'$ and $\mathcal{T}$, so: $$\delta[n] = h[n] \ast h'[n] = h'[n] \ast (\delta[n] - \delta[n-1]) = h'[n] - h'[n-1]$$
So for $n=0$ I have: $$h'[0] - h'[-1] = 1$$
And for $n \ne 0$ I have: $$h'[n] - h'[n-1] = 0 \Leftrightarrow h'[n] = h'[n-1]$$

Still, this isn't enough to solve for $h'$ if I am not mistaken, as for example the following different signals would satisfy this:

$... 4 \ 4 \ 4 \ 4 \ 4 \ 5 \ 5 \ 5 \ 5 \ 5 \ldots$, Here: $h'[-1] = 4$, $h'[0] = 5$
$... 2 \ 2 \ 2 \ 2 \ 2 \ 3 \ 3 \ 3 \ 3 \ 3 \ldots$, Here: $h'[-1] = 2$, $h'[0] = 3$

Can anyone give me a hint what I am missing here or perhaps point me to an error or in the right direction?


texts: I have a multispectral satellite image and I would like to "transform" each of the image bands into one with integer values. For example to stretch? the image between 0 and 255. Not sure how to proceed with this. I also thought I could make a histogram and replace each pixel value with the bin edge in which it falls into. I'm not sure. I'm currently working in python. Any thoughts on this?


texts: How to estimate filter coefficients of an all-pole system?

texts:
Suppose we have measured the frequency response of a system that is known to be all-pole; measuring impulse response is not possible. What are the methods available, if any, to estimate the coefficients of the underlying IIR filter?

EDIT: Frequency response exists only as measured data, so closed mathematical expression is not available, nor is the phase information.


texts: What Distance Measure Is Proper for Image Patches?

texts:
I have a bunch of image patches extracted from two images.

I want to find a proper measure of similarity between each patch from first images to patches from second image.

I have made vectors from each patch to calculate distance/similarity.
I have tried to use euclidean distance and closest vector by phase, but these two measures do not perform well when reconstructing first image by replacing each patch by the most similar patch of the second image.

Is there any well suited measure of similarity for small patches of two images?


texts: Difference between these two digital down conversion methods

texts:
I can think of two possible methods of digital down conversion. One seems superior to the the other for most cases, but I'd like to get some DSP experts' comments on the practical differences between these two.

For both methods described below, I start with the assumption that the $x(t)$ input is already constrained to the band of interest, to avoid the consideration of antialiasing. Antialiasing filters may be necessary in practical applications.

In method A, the in phase component is generated by mixing $x(t)$ with $\cos(\omega t)$ and low pass filtering out the undesired high frequency images. Since there is excess bandwidth in the signal you can use a decimating filter instead of a straight lowpass, shown in my figure below as a blue box around a low pass and a decimator. The quadrature component is similarly obtained by mixing with $\sin(\omega t)$ and then filtering. As a formality I show multiplying the quadrature component by $j$ and adding it to the in phase component.

In method B, convolution is used to perform the discrete Hilbert transform on $x(t)$. This is the reconstruction of the missing quadrature component in $x(t)$. Downconversion is accomplished by mixing the resultant complex signal with a complex sinusoid, no unwanted images generated! Depending on the bandwidth of the baseband signal, it may be smart to add a decimating filter to reduce the data rate, but the output here is usable without that filter. That is, it is not aliased (provided the input signal is limited to the appropriate band).

The primary constraint that I can see in method A is that you can't downconvert by small amounts relative to your bandwidth. Say, for example, you had a 1kHz bandwidth in your signal of interest. If you want to downconvert by 100Hz, method A will not allow for that. The band of the images generated in mixing will intersect the band of the signal of interest. I've seen vector analyzers with analog front ends that can use double down conversion to reject images. A similar method could be used in DSP of course. The problem isn't unsolvable, it just isn't solved by method A.

Method B appears to work around the problems with method A. The construction of a complex signal from the real signal essentially eliminates all negative frequency components, and mixing with a complex sinusoid ensures no images are created, so the only signal present in the output is the signal desired. No filtering is required. However, the convolution required to do the discrete Hilbert transform could be expensive. I would generally choose method A.

My question is: are there advantages / drawbacks of these methods I haven't thought of? What are some practical considerations I should keep in mind before choosing a method?


texts: Is magnitude information enough to reconstruct an audio signal

texts:
I have used an MMSE STSA estimator to obtain the magnitude of an audio signal. The original signal is combined with white noise and I used an algorithm given in an old research paper by Ephraim and Malah to estimate the amplitude of the clean signal. Now I just have the magnitude information and I don't have the phase. The question is: is this enough to reconstruct an audio signal using ifft. If not can I use any arbitrary phase like the white noise phase? Please let me know

Many thanks in advance


texts: How to get general projection of rect function?

texts:
To be more clear I want to get a general expression for

PƟ(t)= ∫∫rect(x)rect(y)δ(xcosƟ+ysinƟ-t)dxdy

I also want to find particular projection for Ɵ=0 and Ɵ=45 and also Fourier Transform of the projections at 0 and 45. I kinda new at image processing and projections, thus I could not figure out a way.

Thanks in advance.


texts: I have two grayscale images, say A and B, and I want to determine, which one is more contrast. For that I'm calculating an RMS contrast as it is described in Wikipedia. The results are about RMS_A = 0.4 and RMS_B = 0.1. Does it mean that image A is more contrast as its gray pixels are nearer to black/white values? How should I interpreter RMS value to determine contrast?


texts: FFT of High Frequency and Low Frequency?

texts:
When we plot FFT in a matlab it doesn't ask for Frequency , all we need to input the sequence and it simply shows the output using stem command.
If i want to compute two different signal using FFT, one with High and the other with Low frequency what should i suppose to do ?
and what will be the difference between the two of them ?

Here is simple code for the FFT calculation ,

x=[1 2 1 0]
n=0:3
y=fft(x)
subplot(2,2,1)
stem(n,y)


texts: Estimating a discrete summer with constrained input bandwidth

texts:
I have a discrete-time system which can be described as:

$$
Y_m = \sum_{r=-N_g}^{R-1+N_g} c_r x[R(m-1) + r]
$$

The unknowns are $c_k$ but I know that they have the following approximate behavior:

$$
c_r \approx \cases{1 \quad 0 \leq r < R \\ 0 \quad \mathrm{else}}
$$

From this, one may recognize that this is nothing more than a discrete summer which obtains the sum of $R$ Nyquist rate samples of $x[n]$. However, in contrast to a normal FIR filter, the output sums are given at a rate $R$ times lower that the one of $x[n]$.

Now my issue is that I do not know $c_k$ exactly. For that reason I want to estimate them. Rewriting in vector matrix form, the Least Squares estimate looks like:

$$
\mathbf{Y}=\mathbf{X} \mathbf{c} \\
\mathbf{c} = (\mathbf{X}^{\#} \mathbf{X})^{-1} \mathbf{X}^{\#} \mathbf{Y}
$$

The success of this method boils down to the conditioning of the matrix $\mathbf{X}$. If $x[n]$ is white, the matrix $\mathbf{X}$ has a low condition number and a rank that is equal to the number of unknowns: $\operatorname{rank}(\mathbf{X}) = R+2N_g$. This is called driving the unknown system with a persistently exciting input. Using simple MATLAB, I can see that the bandwidth of $x[n]$ must be at least $f_s/2 \cdot 0.9$ where $f_s$ is the Nyquist rate of the sequence $x[n]$. Under this condition, the estimation works as expected.

Due to physical constraints (ultimately, this is just a model of a physical system!), I cannot make $x[n]$ white, it needs to be bandlimited to $< f_s/2 \cdot 0.8$. But with such a signal, the rank of the matrix is much lower than the unknowns.

Is there any way to estimate $c_k$ with a non-white $x[n]$ sequence?
Potentially, only for a certain set of input signals (e.g., signals, that are also bandlimited, ignoring the behavior around $f_s/2$ ?
Can it be useful that I know the basic behavior of $c_r$ ?


texts: How to analyse anti-causal discrete transfer function using matlab?

texts:
Consider a discrete transfer function that represents an anti causal filter such as a derivative filter:
$$H(z) = (-z^{-2} -2z^{-1} +2z +z^2) (1/8T)$$

Where T is the sampling period.
Normally in MATLAB, I use fvtool and enter coefficients but these coefficients are relevant to negative exponentials only, i.e. for a causal filter.

I want to view the magnitude response, phase response, pole zero map of an anti-causal filter in MATLAB. How is that possible?


texts: Purple Noise Modeling / Differentiated White Noise (Kalman Filtering)

texts:
I am designing a Kalman Filter for a signal which features a certain kind of noise and I do not know how to model it properly in the filter.
The noise is constructed from a white noise source, called $w(t)$
by taking the difference of the current noise-sample and the last noise sample ($\Delta t$ is the sample-time):

$\eta(t) = w(t) - w(t-\Delta t)$

From my current point of view this is basically high-pass filtered noise or some kind of differentiated white noise (I read purple noise somewhere) with a sinusiodal power spectral density.

Question 1:
Do you agree on my thoughts above?

Question 2:
How can this noise be modeled in terms of a Kalman Filter?
Would a Markov-model be suitable? How to express the color of this noise in a Kalman Filter?

Regards


texts:
I know that group delay is defined as :

$$\begin{align} \tau(\omega) = grd [H(e^{j\omega})] &=-\frac{d[\arg[H(e^{j\omega})]]}{d\omega}\end{align}$$

But if what is available is group delay itself, is it possible to recover the phase?

According to Discrete-Time Signal Processing 3rd edition by Alan Oppenheim equation 5.86 the expression would look like this:

$$\begin{align} \arg[H(e^{j\omega})] = -\int\limits_{0}^{\omega} \end{align} grd [ H(e^{j\phi})]\, d \phi + \arg[H(e^{j0})] $$

and my doubt is:

Where did the second term $$\arg[H(e^{j0}) ]$$ came from?


texts: An effective algorithm for convolution of very many long signals

texts:
I need to calculate the following linear convolution

$$ y[n] = h_1[n] \circledast h_2[n] \circledast h_3[n] \circledast \cdots \circledast h_{k-1}[n] \circledast h_k[n]$$

where $k$ exceeds $5000$ and the length of each $h_i[n]$ exceeds $100,000$.  I tried the method which is based on the circular convolution/DFT but it does not work because each $h_i[n]$ must be zero padded so DFTs then takes too much time and memory. Are there any other algorithm available?


texts: Hints and helps to implement an accurate wavelet packet transfrom

texts:
I want to implement a Time-Freqency analysis by wavelet packets. I don't know what level and sampling rate should I use. I've tried different value for level(3,4,5,6,). but I couldn't get an accurate frequency identification. here is my little scrip and results, hope one can say whats the problem, Thanks in advance.

f=30;N=300;                            %df:final time,N:segments
t=linspace(0,df,N); dt=t(3)-t(2);      % dt: time intervals
fs=ceil(inv(t(2)-t(1)));              % fs: sampling freq
fn=[0.4 0.6 1.0];E=[0.03 0.05 0.1];   %fn:freq, E:damping
 A=[1.0 2.0 3.0];wn=2*pi.*fn;a=E.*wn;wd=wn.*sqrt(1-E.^2);  % A:Amplitude.
y=A(1)*exp(-a(1)*t).*sin(wd(1)*t)+A(2)*exp(-a(2)*t).*...
      sin(wd(2)*t)+A(3)*exp(-a(3)*t).*sin(wd(3)*t);
level = 3;
wpt = wpdec(y,level,'coif5');
[Spec,Time,Freq] = wpspectrum(wpt,fs,'plot');

here are result for Level 3,4,5,6 decompositions.

as you can see in the pictures, It doesn't show fn=[0.4 0.6 1.0] values accurately.
Then I have run another little scrip to see which frequencies,WP identified.

absSpec=abs(Spec);
modI=max(absSpec');          % get  max |Spec| correspond to freq of signal
modI=modI/max(modI);         % normalize each column 0-1
fig=1;fig=fig+1;figure(fig)
Freq=fliplr(Freq);            % the frequency order must be reverse. idk why, i fount-
                              % out by trail and error. very strange but true.
plot(Freq,modI)              % Plot all modes(freq).
grid on

Here is the frequency plot for level 3,4,5,6.

As you can see it doesn't show [0.4,0.6,1]Hz frequencies at level 3,4,5. but in level 6 it shows something(why?). but not accurate though. it represent a range of frequency specially for 0.6 and 1 Hz.(the peak is flat) why this is happening?
moreover, I got the plot for level 8,9,10,11 too. in here I can see a better results, but still it represent 2 or more peaks for each frequencies(why?)(plot is small here so you may barely see it).

so my questions are:

what sampling rate($f_s$)should i select in this case?(now it's
10)
why frequency plot doesn't show any result in level 3,4,5?
It shows some result in level 6 for all frequencies. but what is
its range? I can understand it is representing a range of
frequencies for every components(0.4,0.6,1Hz) of my signal.
Afterward, what is the error for detected range in this level(or any
level)?
In level 8,9,10,11 frequency plot has a narrower range so its
more accurate, but what is other smaller peaks around detected
frequencies?
for a time series with size of M what is WPT computational
complexity in terms of big O(), I thought its faster than CWT but it
is definitely not (I just can feel it by Tic_Toc function). so would
anyone compare WPT with CWT in computational complexity (big O)
terms?

Thanks in Advance for your ideas and answers.


texts: Difference between MATLAB zp2sos and python zpk2sos functions


texts:
I have a filter, with coefficients given by
$$
W=\sin((2*\pi*2*F_f*k/F_s)/(2*\pi*2*F_f*k/F_s))*h[k]
$$, where $F_f$ is the filter reference frequency, $k$ runs from $-N/2$ to $N/2$ ($N$ is the order of the filter), $F_s$ is the sampling frequency and $h[k]$ is a hamming window. For any value of N what will be the attenuation of magnitude encountered by a signal passed through it.


texts: Difference between convolving before/after discretizing LTI systems

texts:
Suppose I have transfer functions for two continuous causal linear-time invariant (LTI) systems:
$F_1(s)$ and $F_2(s)$.

Let $D\left\{\cdot\right\}$ denote the function that maps a transfer function from the continuous-time domain to the discrete-time domain (as another transfer function via the $z$-transform) using ZOH discretization.

What is the relationship between: $$D\left\{F_1(s) F_2(s)\right\} \quad \longleftrightarrow \quad D\left\{F_1(s)\right\}D\left\{F_2(s)\right\}$$

In other words, what is the difference between discretizing the convolution  of two transfer functions versus convolving their discretization?
Some brief experimentation has shown they are "almost" equal, but why?
Can we be precise about this difference?

As an aside,

Could someone point me to a reference on other properties of discretization, e.g., linearity $D\left\{a_1F_1(s) + a_2F_2(s)\right\} = a_1D\left\{F_1(s)\right\} + a_2D\left\{F_2(s)\right\}$?


texts:
the following MATLAB code should plot the iFFT of $-i \pi sgn(f)$ which has to be $\frac{1}{t}$ in continuous time. As you can see, the output kinda looks like a 1/x function, but is modulated.
How can this phenomena be explained and is there a way to get the desired $\frac{1}{t}$ function?

clc;
clear;
close all;

L=1000;
fStart=-100;
fEnd=100;
df=(fEnd-fStart)/L;
f=linspace(fStart,fEnd-df,L);

X=-1i*pi*sign(f);
x=ifft(fftshift(X));

figure()
ax1=subplot(3,1,1);
ax2=subplot(3,1,2);
ax3=subplot(3,1,3);

stem(ax1,imag(X));
plot(ax2,fftshift(real(x)));
plot(ax3,fftshift(imag(x)));

title(ax1,'spectrum');
title(ax2,'Re(iFFT)');
title(ax3,'Im(iFFT)');

Edit: for clarification: I especially wonder, why every other sample in the result is zero:


texts:
Is it possible to blindly estimate a coarse carrier frequency offset of an OFDM signal where $\lvert f_{\rm off}\rvert > \frac{\Delta f}2$, where $\Delta f$ =  subcarrier spacing? Blind is defined as no a-priori knowledge of the OFDM signal parameters. However, if the useful period, Tu, is required then that can be assumed to be known.

I am wondering if there is a technique for coarse frequency estimation similar to the fine frequency estimation outlined by van de Beek et al $^1$.

$^1$J. J. van de Beek, M. Sandell, and P. O. Borjesson, “ML estimation of time and frequency offset in OFDM systems,” IEEE Trans. Signal Process., vol. 45, no. 7, pp. 1800–1805, Jul. 1997.


texts: Detecting the original bandwidth of upsampled recordings

texts:
I have a set of recordings sampled at 16khz, but some of them are upsampled from 11khz and some from 8khz. Quite some of the recordings are noisy. What is the best way to classify the recordings according to the original bandwidth?

I know some people are using GMM classifiers for detection. A paper on the subject is DETECTING BANDLIMITED AUDIO IN BROADCAST TELEVISION SHOWS

I've come to the following algorithm: calculate energy for 40 frequency bins for the full range from 0 to 8khz for 10ms frames, then normalize by maximum energy, then calculate maximum of normalized energy over sequence of frames. Bandwidth is detected by comparing maximum normalized energy in high-order bins compared to maximum normalized energy in middle bins.

Any better ideas?


texts: Decision directed TED (timing error detector). What does it mean?

texts: I am reading about symbol rate recovery and there is a term that I didn't grasp - 'decision directed' TED. Can somebody explain what does it mean or provide me with accurate reading? Thanks


texts: I want a relatively easy to implement but relatively reliable algorithm for detecting beat locations in an audio file.What I am aiming to do is get the location of the closest beat to a particular time instant i.e if given current playtime, how can I get the location of next beat.Does anyone know of any algorithms that could help me solve this problem?


texts: How to reduce latency under mean filter for high noise?

texts:
We have a position sensor that under some conditions receives some high frequency noise.  We can eliminate that very well with a simple mean filtering.

Unfortunately this causes too much lag when there is not any noise (at the beginning of the graph at s=90 to 95). Thus we have been experimenting with Savitzky-Golay filtering:

However, the results are not that great. The filtering under heavy noise is worse than mean filtering (Or the latency reduction is not that good). What else can we do?
Some more background:

The problem is an online problem, so we cannot retrospectively analyze the data and for instance shift the mean by half a period.
The 'regular' movement is performed by a human moving an object, the 'noise' is from a high frequency vibration.


texts: Transfer function of double exponential smoothing?

texts:
Double exponential smoothing
a.k.a.
Holt-Winters smoothing
tracks level and trend of a time series in coupled IIRs:
$\qquad$ In:  $Y_t$, t = 0 1 2 ...
$\qquad$ State: $L_t, T_t \quad$ -- level and trend
$\qquad$ Out: estimate $\hat{Y}_{t+1}$
$\qquad$ Parameters: a, b  (a.k.a. $\alpha, \beta$)
$\qquad$ Step equations:
$\qquad$ $\qquad L_t = a Y_t + (1 - a) (L_{t-1} + T_{t-1})  \quad$ -- level
$\qquad$ $\qquad T_t = b (L_t - L_{t-1}) + (1 - b) T_{t-1}  \quad$ -- trend
$\qquad$ $\qquad \hat{Y}_{t+1} = L_t + T_t$

How can I calculate the transfer function of this smoother, for given $a$ and $b$ ?

(Two possible approaches:

manipulate the step equations into ratio-of-polynomials form
brute force iterate with input $e^{2 \pi i f t}$: converges slowly for small $a, b$.

The first of these is beyond me, the second unsatisfactory.)


texts: Having converted a grayscale image to binary using `im2bw()`, is the converse operation possible?

texts: I have converted a gray-scale image to binary using im2bw() now i wanna do the opposite. How can I do that in MATLAB?


texts: How to design an inverse for a Butterworth low-pass filter?

texts:
I have a signal that has a spectral shape very similar to the transfer function of a second-order lowpass filter (butterworth design).

I would like to whiten the spectrum of this spectrum in the time domain. For that, I need a filter that reverses the effect of the second-order lowpass filter.

I can calculate the inverse filter of the second-order lowpass by exchanging the A- and B-coefficients, but this will amplify high frequency components to infinity.

How do I design a filter that has the inverse transfer function of a second-order butterworth design but does not amplify to infinity at $f_s/2$? The transfer function should have the same corner frequency, but smoothly rise with 12 dB per octave.

Is this even possible? Any pointers as to how to design such a filter would be helpful!


texts:
I've tried to implement STFT in functional way in Python. I've got code such as below:

def genSpectrums_dft(x, w, N, dft_func = fp.fft):
    hN = N / 2
    M = len(x)
    hM1 = int(math.floor((M + 1) / 2))
    hM2 = int(math.floor(M / 2))
    wM = w(M)
    w_sum = sum(wM)
    xw = x * wM / w_sum
    dft_buffer = np.concatenate(
        (xw[hM2:], np.zeros(N - hM1 - hM2), xw[:hM2])
    )
    X = dft_func(dft_buffer)
    mX = 20 * np.log10(abs(X[:hN]))
    pX = np.unwrap(np.angle(X[:hN]))
    return mX, pX

def genSpectrums_stft(x, w, M, N, H, dft_func = fp.fft):
    hM1 = int(math.floor((M + 1) / 2))
    hM2 = int(math.floor(M / 2))
    x_ext = np.concatenate((np.zeros(hM2), x, np.zeros(hM2)))
    lv = np.arange(len(x) / H + 1)
    wM = w(M)
    w_sum = sum(wM)
    xlv = np.array([x_ext[l * H:l * H + M] / w_sum for l in lv])
    return tuple(
        np.transpose(np.array(X))
        for X in zip(
            *[genSpectrums_dft(xl, w, N, dft_func) for xl in xlv]
        )
    )

Full working code is available here (Python, numpy, scipy and matplotlib required to run it).

When I run my code for any wave, STFT magnitude spectrum seems to be quite good, but STFT phase spectrum always looks like this:

But as far as I know it should rather look like this:

Both these spectrums have been generated from the same wave and the second parameters. The second one has less frequencies. Why my phase spectrum has these strange vertical lines?

I've got a few remarks rewarding my code:

I am almost positive that finding dft spectrums works correctly
Though the spectrums distinct, when I inverse the process, I get something really similar to the input

If you are not familiar with functional programming, I don't mind if you correct my mistakes in structural way, but the most important thing for me is to understand the problem with my code.


texts: Why are edges in spatial images represented as edges in their Fourier transform image?

texts:
Here is a well-known image and its Fourier Transform (magnitude).

If I understand correctly the theory behind the FFT, each pixel in the FFT image represents a certain 2D sine wave with frequency depending on the distance from the center of the image, and orientation depending on angle with the horizontal. Intensity of these pixels indicates the coefficients with which each sine wave is added, which, combined with (hidden here) phase information, gives the original image back if we do an inverse Fourier transform.

While I fail to grasp is that when treating real images such as the ones presented, why are edges in the images visible as edges in the frequential domain ?
In this example, there is a diagonal line in the FFT image (let's forget about the vertical and horizontal line which I think are artefacts based on the way FFT is computed, needs a periodical image, etc). This diagonal line is probably caused the girl's hat. But as I understand it, the line in FFT domain means a sum of sines oriented in the same way but with different frequencies. How does that result in an edge when we convert back ? Since edges are high frequency information, wouldn't an edge be represented by one very bright point in the FFT ?
Does it have to do with adding different sines so that they cancel each other out on some portions of the image ? Does the phase image has anything to do with it ?

A more compelling example might be the following image set :


texts: How to define SNR for multiple signals in noise?

texts:
I have to estimate the number of signals present in a measurement contaminated by additive noise given $n$-dimensional snapshot vectors $\bf x$, modeled as $ \bf x = \bf A \bf s + \bf z $ where $\bf s$ is a $k \times 1$ vector representing $k$ different signals. $\bf A$ is an $n \times k$ non-random matrix and $\bf z$ is an $n \times 1$ noise vector. (This model is common in array processing  problems.)

How should I define the SNR—per signal, averaged or otherwise?


texts: How do I resample an audio signal without altering its pitch?

texts:
I'm working on concatenating videos again, this time with different playback speeds. The video stream part is easy enough, but the audio part is what's stumping me at the moment.

Currently, I've managed to do it by just multiplying the timestamp by a factor inversely derived from the desired playback speed, so that if the playback is to be half speed, the timestamps are multiplied by two.

Not only that, I've also altered the audio samples (in the form of a byte array) so that when the playback is half speed, every sample (of four bytes) is repeated, and when the playback is double speed, every second sample is discarded.

My problem is that this method also alters the pitch of the resulting audio track - the half-speed playback lowers the pitch, the double-speed playback raises it.

How can I resample these samples without altering their pitch?


texts: `SSIM(YUV)` vs. `SSIM(Y)`: Values Are Very Different

texts:
I'm seeing some weird results when trying to compare image qualities in a few different color spaces. Have a look at my graph:

I want to compare with another publication, which reported PSNR(YUV) and SSIM(Y). However, I don't understand how I can have SSIM(YUV)=0.9 for the highest quality point on the right, but only have SSIM(Y)=0.7 for the exact same image. How can there be such a big difference? Is this normal? Or am I doing something wrong?

For completeness, here are reference image and the image under test:

And the MATLAB code I use for converting RGB to YUV:

function YUV=rgb2yuv(RGB)
R = RGB(:,:,1);
G = RGB(:,:,2);
B = RGB(:,:,3);

Y = 0.299   * R + 0.587   * G + 0.114 * B;
U = 0.5 - 0.168736 * R - 0.331264 * G + 0.5 * B;
V = 0.5 + 0.5 * R - 0.418688 * G - 0.081312 * B;
%   ^^^ I tried to use 0.0 as well, but SSIM values are identical

YUV=cat(3,Y,U,V);
end


texts:
I'm attempting to code Pisarenko's method for some time-series data that I have.

I've had a look at a couple of papers, which explain in detail the theory behind it, but I'm having some trouble figuring out how to convert the math into an algorithm for either C or Python.

Is there any information on how one can do that? A pseudocode algorithm, or some hints on which direction to start etc?


texts: Extracting features from a signal: a 21st century approach?

texts:
I have a discrete signal of random length (let's say between 500 and 800 in time). This signal has an important information which I want to collect and to put it together in a feature vector to train a classifier. The signal also comes in random phase and in random scale (it depends on the distance to the sensor).

So there are some constraints: the feature vector has to be invariant to scale and phase, it has to be a sparse representation of the information contained in the signal and it has to deal with the problem of different signal length.

I have followed a classic approach to solve this problem:

I calculate the Fourier coefficients of the signal, a[n] and b[n]. I take as feature vectors the most significant coefficients for my problem (I fix this parameter experimentally). To solve the problem of scale and phase invariance I consider only the normalized magnitude r[n] = s[n]/s[0] where s[n] = sqrt(a[n]^2+b[n]^2).

This solution actually works quite well, but I am sure that there has to be something more "up-to-date" to solve that. I would really appreciate any alternative approaches to solve this problem or some references to read.

Thank you.


texts: Impulse and step response of a system that contains cascaded second order sections

texts:
I have to design an 8th-order IIR filter with the following transfer function magnitude response as a specification:

I use "ellip" function from matlab to calculate the unquantised coefficient set:

Then i implement the filter using a cascade of second-order sections using "tf2sos" matlab function:

And then i quantize these coefficients:

After all these steps i take the two matrices below:

I choose the filter representation for every second-order instance to be like below (Direct Form 2):

circle with "X" is a multiplication
circle with "Σ" is an addition
and the triangle is a delay element

After that i have done all of these i write some matlab code to model the structure above:

%Implementation of a second-order filter structure
function[y,state,pointer]=Second_Order(x,a,b,state,pointer)

%The filter routine when fed with an input sample produce one output
%sample each time it is called

%x=input value
%a=array containing feedback coefficients[1,a1,a2]
%b=array containing feed forward coefficients[b0,b1,b2]
%state=array(of length 2) containing old input values
%y=output value

%compute output value
y=b(1)*x+state(pointer+1);
pointer=rem(pointer+1,2); %update pointer in modulo form

%Update states
state(pointer+1)=b(2)*x-a(2)*y+state(pointer+1);
pointer=rem(pointer+1,2);

%Overwrite oldest sum with b(N-1).x
state(pointer+1)=b(3)*x-a(3)*y;
pointer=rem(pointer+1,2);%Increment pointer modulo-(N-1)

At the end, and the point that i don't know how to do it properly is how i can display the impulse and step response of the whole filter and not for every second order section.

For example if i write something like that, i think i take the impulse response of the first second order section. How i can take the impulse resonse of the whole system?

%Initialization
[b,a]=ellip(4,1.3,60,2*[0.2,0.3]);
[sos,q]=tf2sos(b,a);
sosq=Quantize('round',sos,2^6);

bf=sosq(1,1:3);
af=sosq(1,4:6);

%impulse invariant design
impsv=zeros([1,2]);%impulse response state variable
stpsv=zeros([1,2]);%step response state variable
len=20; % run length
impy=zeros([1,len]);%impulse response result
stpy=zeros([1,len]);% step response result
impp=0;
stpp=0;

% Simulation of filter operation

x=1;
for k=1:len
  [impy(k),impsv,impp]=Second_Order(x,af,bf,impsv,impp);
  [stpy(k),stpsv,stpp]=Second_Order(1,af,bf,stpsv,stpp);
  x=0;
end


texts: Can't seem to get welch function in scipy work correctly

texts:
I'm trying to learn DSP through using Python's scipy package.
I have measure some 200Hz signal from a machine. Now I want to inspect the spectral density of the signal. Here's the signal plotted out:

As you can see, the signal is obviously not a low frequency signal. However, I still get the PSD graph like below.
Here's the code I have been using:

def plot_and_save(sample,i=""):
filepath = sample + "/" + "merged_data.csv"
df = pd.read_csv(filepath)

sns.set_style("darkgrid")
sns.set_colodef plot_and_save(sample,i=""):r_codes("dark")
fig1, axes = plt.subplots(3,2)
fig1.figsize = [1920,1080]

axes[0,0].plot(df["Time"], df["Voltage"], marker="",color="blue")
axes[0,0].set_title("Plot of Voltage against time")
axes[0,0].set_xlabel("Time (us)", ha= "right", x=1.0)
axes[0,0].set_ylabel("Voltage (V)")

axes[0,1].plot(df["Time"], df["Current"], marker="",color="red")
axes[0,1].set_title("Plot of Current against time")
axes[0,1].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[0,1].set_ylabel("Current (A)")

sampling_freq = 1e6
vol_f, vol_psd = welch(df["Voltage"].to_numpy(), fs=sampling_freq)
axes[1,0].plot(vol_f, vol_psd, marker="",color="green")
axes[1,0].set_title("Plot of Voltage's Power Spectral Density")
axes[1,0].set_xlabel("Frequency", ha= "right",x=1.0)
axfilepath = sample + "/" + "merged_data.csv"
df = pd.read_csv(filepath)

sns.set_style("darkgrid")
sns.set_color_codes("dark")
fig1, axes = plt.subplots(3,2)
fig1.figsize = [1920,1080]

axes[0,0].plot(df["Time"], df["Voltage"], marker="",color="blue")
axes[0,0].set_title("Plot of Voltage against time")
axes[0,0].set_xlabel("Time (us)", ha= "right", x=1.0)
axes[0,0].set_ylabel("Voltage (V)")

axes[0,1].plot(df["Time"], df["Current"], marker="",color="red")
axes[0,1].set_title("Plot of Current against time")
axes[0,1].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[0,1].set_ylabel("Current (A)")

sampling_freq = 1e6
vol_f, vol_psd = welch(df["Voltage"].to_numpy(), fs=sampling_freq)
axes[1,0].plot(vol_f, vol_psd, marker="",color="green")
axes[1,0].set_title("Plot of Voltage's Power Spectral Density")
axes[1,0].set_xlabel("Frequency", ha= "right",x=1.0)
axes[1,0].set_ylabel("PSD")

cur_f, cur_psd = welch(df["Current"], fs=sampling_freq)
axes[1,1].plot(cur_f, cur_psd, marker="",color="yellow")
axes[1,1].set_title("Plot of Current's Power Spectral Density")
axes[1,1].set_xlabel("Frequency", ha= "right",x=1.0)
axes[1,1].set_ylabel("PSD")

axes[2,0].plot(df["Time"], df["Current"]*df["Voltage"], marker="",color="red")def plot_and_save(sample,i=""):
axes[2,0].set_title("Plot of power against time")
axes[2,0].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[2,0].set_ylabel("Watt")

sns.scatterplot(x="Current", y="Voltage", data=df, ax=axes[2,1], color = 'r', edgecolor="None", linewidth=0)
axes[2,1].set_title("Plot of Voltage against Current")
axes[2,1].set_xlabel("Current (A)", ha="right", x=1.0)
axes[2,1].set_ylabel("Voltage (V)")

fig1.suptitle("Sample number: " + str(i+5),fontsize=16)
# plt.savefig("formatted/"+date+"/figure/"+str(i+5)def plot_and_save(sample,i=""):".svg",format="svg")
plt.show()es[1,0].set_ylabel("PSD")

cur_f, cur_psd = welch(df["Current"], fs=sampling_freq)
axes[1,1].plot(cur_f, cur_psd, marker="",color="yellow")
axes[1,1].set_title("Plot of Current's Power Spectral Density")
axes[1,1].set_xlabel("Frequency", ha= "right",x=1.0)
axes[1,1].set_ylabel("PSD")

axes[2,0].plot(df["Time"], df["Current"]*df["Voltage"], marker="",color="red")
axes[2,0].set_title("Plot of power against time")
axes[2,0].set_xlabel("Time (us)", ha= "right",x=1.0)
axes[2,0].set_ylabel("Watt")

sns.scatterplot(x="Current", y="Voltage", data=df, ax=axes[2,1], color = 'r', edgecolor="None", linewidth=0)
axes[2,1].set_title("Plot of Voltage against Current")
axes[2,1].set_xlabel("Current (A)", ha="right", x=1.0)
axes[2,1].set_ylabel("Voltage (V)")

fig1.suptitle("Sample number: " + str(i+5),fontsize=16)
# plt.savefig("formatted/"+date+"/figure/"+str(i+5)".svg",format="svg")
plt.show()

Thank you so much for the help!


texts: Can't get consistent magnitude response in MATLAB

texts:
I have designed a notch filter with MATLAB filter designer app. Then I generated the filter design function which gives me the filter object:

function Hd = filter_notch

Fs=500; Fnotch=50; BW=1; Apass=1;

[b, a] = iirnotch(Fnotch/(Fs/2), BW/(Fs/2), Apass);
Hd     = dfilt.df2(b, a);

Now, I use freqz() to plot the response. However I get 2 different responses.

Mehtod 1. Direct plot with freqz(filter_object)

>> freqz(filter_notch);

Method 2. Manual plot with h=freqz(filter_object, Physical_frequencies, F_sample)

>> fs=500;
>> f_range=linspace(0,fs/2,2048);
>> h=freqz(filter_notch,f_range,fs);
>> hdb=20*log10(abs(h));
>> plot(f_range,hdb)

you can see the notch attenuations are different with above 2 methods, roughly -25dB vs -15dB. Why is that?

From MATLAB documents:

 h = freqz(___,f,fs) returns the frequency response vector, h, at the
 physical frequencies supplied in f.


texts:
Suppose $g[n]$ is periodic with fundamental period $N$ and $f[k]$ being its Fourier coefficients i.e. $$ f[k] = \frac{1}{N}\sum_{n=<N>}g[n]e^{-j\frac{2\pi}{N}nk}$$ In more convenient notation $$g[n]\leftrightarrow f[k]$$  It's easy to prove $$f[n]\leftrightarrow \frac{1}{N}g[-k]$$ This is called "Duality in the discrete-time Fourier series". My problem is applying duality in the case of $$x[n-n_0] \leftrightarrow a_ke^{-jk(\frac{2\pi}{N})n_0}$$ And $$\sum_{r=<N>}x[r]y[n-r]\leftrightarrow Na_kb_k$$
I think for the first case duality should gives$$x[n]e^{jk(\frac{2\pi}{N})m} \leftrightarrow \frac{1}{N}a_{m-k} \tag{1}\label{*}$$ and the second gives $$x[n]y[n] \leftrightarrow \frac{1}{N^2} \sum_{l = <N>}a_lb_{-k-l} \tag{2}\label{**}$$
Edit:
Let $g[n] = \sum_{r=<N>}x[r]y[n-r]$ and $f[k] = Na_kb_k$. So according to the mentioned property $$Nx[n]y[n] \leftrightarrow \frac{1}{N}\sum_{r=<N>}a_rb_{-k-r}$$
Dividing by $N$ leads to $\eqref{**}$.
Similarly if $g[n] = x[n-n_0]$ and $f[k] = a_ke^{-jk(\frac{2\pi}{N})n_0}$ I think then $$x[n]e^{-jn(\frac{2\pi}{N})n_0}\leftrightarrow a_{-k - n_0}$$
Let $m = -n_0$ and we have \eqref{*}.
Obviously these results are wrong but I don't know what's my mistake.
This is taken from Signals and systems by Alan V. Oppenheim:


texts: Carrier Tracking Loop: Tracking noise vs Phase Noise Trade-off

texts:
I would like to compare two different implementations of carrier tracking loop which differ in the way phase detection happens in the PLL. I am using PSK modulation for the data.
In the first implementation, a mixer is used as a phase detector and the implementation looks like this:

Here, to ensure that the PLL does not track the PSK modulated data, the bandwidth needs to be sufficiently below the symbol rate to ensure that the tracking noise is not too high. This I can understand because if only data (without any noise) is passed through the PLL, low frequency components of data will be tracked and will contribute to noise. This can be seen in the figure below, where the output data from the PLL is high pass filtered because of the Loop:

Another implementation where the phase detector is implemented as a vector multiplier is explained here. In this implementation, I do not clearly see if there would be such kind of tracking noise vs phase noise trade-off. Because if I pass data without noise, than phase detector output will be zero and there would be no filtering of the data by the PLL. The plot for the PLL output with vector multiplier phase detector is shown below and I see that the tracking noise is non-existent even though I chose bandwidth very close to the data rate.
Why would such a trade-off exist for this kind of carrier tracking implementation?

EDIT
For the mixer implementation I am using ideal phase detector:

delta_phi = Kd*(phi_ref - phi_vco) %% Kd is the PD gain

I use complex exponential to implement vector multiplier:

phi_rotation = ref.*vco';  %%% Both and vco are complex exponentials
symbol_estimate = conj(pskDemodulator(phi_rotation));
delta_phi = angle(phi_rotation.*symbol_estimate);  %%% imag would give the same result.

Block diagram Approach 2

Block Diagram Approach 1


texts: What is the acceptable range for bandwidth utilization factor in OFDM systems?

texts: The bandwidth utilization factor is defined by (alpha=bite rate/bandwidth), for OFDM system, what is the acceptable range for alpha?


texts: How to recover $f(t)$ from Fourier Transform of its absolute value $\mathcal{F}|f(t)|$?

texts:
Let the Fourier Transform of a real signal, $f(t)$, be $\mathcal{F}(\omega)$. And the FT of the absolute value of the same signal, $|f(t)|$, be $\mathcal{F}(u)$.

Can $\mathcal{F}(w)$ be recovered from $\mathcal{F}(u)$?

For instance, the FT of $a \cdot \cos(ft)$ returns a spectrum in which the frequency $f$ has amplitude $a$.

Can $f$ and $a$ be recovered from the FT of $a \cdot \cos(ft)$?


texts: How we calculate Precision-Recall Curve?

texts:
As far as I know, precision and recall are two single values. How we can plot a curve from these two single values? I think I should calculate a set of values for each of them but how?

afterwards, the curve can be depicted by using the fact that they are inversely related.

$$
\text{Precision}=\frac{tp}{tp+fp}
$$
$$
\text{Recall}=\frac{tp}{tp+fn}
$$
Where $tp = \text{True Positives}$, $fp = \text{False Positives}$ and $fn = \text{False Negatives}$.

Anybody can explain how we can plot a curve from two single values?

One answer can be found here but I do not believe that I have caught the point well.


texts: Using Scipy.correlate2d output matrix to determine best image match

texts:
I'm struggling to wrap my head around how I can use the output of a correlation matrix to determine if my image is present.
If I am trying to match characters, how can I compare the outputs of these correlation matrix's to determine the best match?
Lets say I have the MNIST database:
if I compare one of the test digits to my labeled training set using correlate2d, how do I know which matrix is the best choice to label my test digit?
for i,j in enumerate(training_digits):

  target_digit = test_digits[1]
  template = training_digits[i]
  corr = signal.correlate2d(target_digit, template, boundary='symm', mode='full')


texts: Implementing real-time signal processing filters

texts: I’m currently working on my capstone for school, and I’m looking for some guidance/tips when it comes to real-time filtering. The nature of the project is to monitor breathing patterns, send the aquired signals wirelessly, and then process/filter the data using a mobile app. My part of the project is the latter, and I would like some tips in terms of coding language/various articles that could guide me in the right direction! Thanks


texts: Effect of sampling frequency to the center of a bandbass filter

texts:
Is it true that if you have some sort of filter, e.g. a digital bandpass filter around the frequency 1000 Hz that is given in the form of a transfer function $H(z)$ and you then change the sampling frequency of the input that the same filter can suddenly have its "center" (of the bandpass) at another frequency?

Which keyword can I look up to learn more about this concept?


texts:
I have made some segments of the image and found the area of interest.I used Sobel vertical detector and morphological operations.

Now I want to segment all number plate characters from the image which has size [51 165]......How would I do it now?Even segmentation from the whole image was a bit challenging task since my technique does not work if the car(since these are number plates) is too near or far or illumination is poor if somebody has suggestions then please guide me on that.

Sample Image from the dataset


texts: fourier transform of smart-phone accelerometer in matlab

texts: I'm new to matlab. I want to process my smart-phone accelerometer data in matlab. I know Matlab let's you connect your phone via USB cable to see accelerometer data in realtime. But according to some reasons I don't want to use this method. Instead I've developed an application that saves accelerometer data. I want this application to produce a ".m" file and use this file to represent x-y-z signals. The problem is that the accelerometer data rate is not constant. for example: "data-1" corresponds to time "0.1", "data-2" corresponds to time "0.2", and "data-3" corresponds to time "0.8". You can see the time interval can change. I need to get fourier transform of these signals, so I can not use the simple matlab vector to represent them. How can I represent these signals? How can I get the correct fourier transform of these signals?


texts: What do the L, M and S pixel values in the LMS color space stand for exactly?

texts:
From what I read in the literature, the LMS color space is an attempt to quantify how the human retina's cone cells respond to light. More specifically, L has to deal with long wavelengths, M with medium wavelengths and S with short wavelengths. But so far I haven't been able to fully grasp the "scientific" meaning of LMS values.

Let's talk about L for example, here's how I interpret its value: A pixel p on the image is just a projection of a point P of the 3D world on the image plane. That world point P sends light radiations of various wavelengths. The L value of pixel p represents how many photons stemming from light radiations of long wavelengths were absorbed by the retina's cones.

To what extent is my understanding correct?


texts: Audio experiment with intermodulation by multiplying two analytic signals

texts:
1) First, I describe intended processing with a basic case:

Let's say on left channel I have a sin wave with frequency F1, and constant amplitude A1.
Then on right channel a sin wave with frequency F2 and constant amplitude A2.
Then I would generate a third sin wave, with frequency F3, and constant amplitude A3. Such third sin wave verify: either  A1≤A3<=A2 or A2≤A3<=A1. Also F3= F1+F2.

I chose to fix A3 as  A3= Sqrt(A1*A2).

The third sin wave would be blended to a small fixed amount  to left
   and right channels (let's say at -30 db) . This would add some intermodulation.

In our particular basic case, I can get the target third sin wave by multiplying analytic signals of left and right. (Context:  Analytic signals are related to Hilbert transform . I use a FIR).

Let's say on left I have  re1 + i*im1, then on right re2 + i*im2 .

Then I can get my third sin wave by computing (two square roots applied on denominator) :

  (re1*re2 - im1*im2) / √√(asq1*asq2)

Asq1 is the "absolute square" for left signal i.e :

  asq1= re12 + im12

Same story with asq2.

You might ask , why not do a straightforward "ring modulation", by multiplying directly both channels, i.e re1*re2 . Because it will introduce an additional harmonic with frequency Abs(F1-F2) that I don't want.

2) Pushing the experiment to true audio content:

Ok, I understand that the audio effect I imagined might be flawed, or that I might be using the Hilbert Transform in a wrong way.

Nevertheless here are few improvements I've made to the implementation of this process:

The frequency response from the Hilbert FIR is not completely flat, and is like a wide bandpass. So I ensured to reproduce same frequency response for the real part of each analytical signal.
I avoid aliasing by upsampling 2x, analytic signals, and then downsampling after multiplication.
The division by the term √√(asq1*asq2) introduces lot of unwanted frequencies. I've mitigated this by applying a gaussian blur to asq0, and asq1.

The division by the term √√(asq1*asq2), remain problematic even after applying some blurring. But I want the amount of distortion to be independant of the level of original signal, so I  keep it.  In other  words, if I adjust the volume, the added distortion will vary in a linear way. Not instantly, but on "average" with time .

The subjective quality of the result, vary with audio content,  the blurring period of asq0 & asq1, and also the length of hilbert FIR. It mostly affect the perception of transients (I remind that the "intermodulation" is blended to original audio, at low level).

Anyways it's an audio experiment.  I could have played instead by adding some second harmonic distortion. It can be pleasing, and it's often used by musicians. But I imagined something different. Hopefully, this looks enough simple.

Now, if you  could point out what's wrong in this processing effect, or how it could be improved, that would be great.


texts: State of the art for automatic rhythm analysis

texts: What are the currently most used and less error prone algorithms in the field of music information retrieval for beat detection and tracking?


texts:
In Hamming's book, The Art of Doing Science and Engineering, he relates the following story:

  A group at Naval Postgraduate School was modulating a very high
  frequency signal down to where they could afford to sample, according
  to the sampling theorem as they understood it.  But I realized if they
  cleverly sampled the high frequency then the sampling act itself would
  modulate (alias) it down.  After some days of argument, they removed
  the rack of frequency lowering equipment, and the rest of the
  equipment ran better!

Are there any other ways to use aliasing as a primary technique for processing a signal, as opposed to a side-effect to be avoided?


texts: Estimation of Input Signal to Obtain the Desired Output Signal for an Unknown Filter

texts:
Suppose $h(n)$ is a finite impulse response which is unknown. We can feed any input signal $x(n)$ into the system and observe the corresponding output signal $y(n)$.

From this, is it possible to estimate the input signal $x^*(n)$ which when passed through the filter gives the desired output signal $y^*(n)$ ?

I think a Least Mean Squares filter could be used here. Could anyone point out how one proceeds doing this ?

Thanks !


texts:
I am doing a project where I need to take the FFT and IFFT of a photo in MATLAB. There is a principle in optics where the intensity distribution from a lens is equal to the Fourier transform of the aperture that the light initially goes through. I'm having issues finding the FFT and IFFT of an ideal source.

A circular aperture should have an intensity pattern called Airy disk, which should result from the Fourier Transform of a circle. Here is a photo of the Airy disk that I'm using in my code:

Taking the inverse Fourier transform of the Airy disk should result in an image of a circular aperture, but all I'm seeing is black when I convert to uint8. If I leave it as type double then I see things that also look like Airy disks in the corners. Here is the code I'm using:

pic=imread('airy.jpg');
pic=rgb2gray(pic); %Convert to grayscale
imshow(pic)
title('Initial');

pic=double(pic); %Convert to double before transforms
pic=ifft2(pic); %Take IFFT of Airy disk
pic=real(pic); %Take real part

figure()
pic=uint8(pic); %If commented out, black screen has Airy disks in corners
imshow(pic);

Does anyone have any ideas on why the IFFT of the Airy disk doesn't return a circle? Thank you


texts: FFT of color images incorporated into an Object Recognition method?

texts:
Is there any instance where Fourier transforms of color images have been used in conjunction with other object recognition method?

Any instance of usage of Fourier transforms in color images? I think there is great potential in the usage of Fourier transforms of color images.


texts: find cross correlation lag by optimization

texts:
I have data from many closely located sensors (geophones). These are often contaminated with coherent noise which I need to remove. Some of the most effective de-noising methodologies work if the noise signals are aligned. To do that I do cross correlation on neighbouring sensor data and find the delay and shift my data accordingly. The problem is the noise source is moving, but slowly. So I would need to recalculate the delay for every new data set. Knowing the initial delay (estimated from the first data set), is there a way, using some optimization methods to find the subsequent delays, i.e. at each step use the previous delay estimate as input (together with new data) to some optimization function and find the delay for the new data set (in Matlab) ?

Thanks in advance


texts:
I have computed 12d-MFCC features on 3 consecutive audio frames, say:

For 1st frame: a1 a2 a3 ... a12
For 2nd frame: b1 b2 b3 ... b12
For 3rd frame: c1 c2 c3 ... c12

How to organize these MFCC features, in a1 a2 ... a12 b1 b2 ... b12 c1 c2 ... c12 or a1 b1 c1 ... a12 b12 c12? Is there any difference between these two organization styles (for audio processing)?


texts: How to transform data to uniform distribution (uniform percentiles)?

texts:
Given the list of data points (normalized in [0,1] range), I plot the histogram of values and compute percentiles (shown as x ticks).

How to find a transformation of data values so the histogram is approximately uniform. Which would, in turn, make percentile values also uniformly distributed.


texts: What causes this type of distortion in a constellation I/Q plot?

texts:
I'm transmitting QPSK through a digital pulse shaping filter (RRC, 35% excess BW, 2 samples per symbol, 64 taps) at a rate of 56Msym/sec.
Receiving it on standard lab equipment (Vector signal transceiver), where it is demodulated with RRC MF and plotting the constellation. At high symbol rates, the phenomenon below starts appearing in the constellation where the points get essentially binned into 4 different locations (looks like some exotic modulation you might call qpsk^2)
One thing of note, this effect is easily taken care of by the equalizer on the RX and reduces the constellation to a much more ideal, low EVM one. Given this, I suspect the effect is related to ISI. One possible thought I had was for example, on the I channel each symbol is perturbed by a neighbor (which could be a + or - symbol) so the effect is two offset bins. Same for Q channel, resulting in the 4 shown. Not sure why this contribution would only come from a single symbol as opposed to many though - in which case I would expect more of a cloud/awgn like effect from adding up contributions of many symbols)
For reference the basic D/A chain looks like:
Generate QPSK IQ -> Filter with RRC -> Halfband Upsamplers -> DAC -> Analog Reconstruction


texts:
I just read a paper yesterday and there is confusion:

I want to know what is reconstruction fidelity term and the theorem of changing the formula from upper to lower one. Is this theorm from functional analysis?

Could anyone help me ?

EDIT: The link of the paper.


texts: Evaluate Fourier coefficients at arbitrary point using Python

texts:
Lets say I have a sinusoidal function $s$ that looks like

import numpy as np
import matplotlib.pyplot as plt

s = lambda t: np.sin(2 * np.pi * t / 13)
s_sampled = np.array([p1(i) for i in range(50)])
plt.plot(s_sampled)

Now I want to take the Discrete Fourier transform of the signal and reconstruct the whole original signal. As the signal only has a period of 13. It should be enough to Fourier transform the first 13 points.

fourier_coefficients = np.fft.fft(s_sampled[:13])

However, using the numpy's inverse transform I can only get back 13 points.

plt.plot(np.fft.ifft(fourier_coefficients))

How can I get the rest of the (now just repeating) signal back?
I have tried to build a function myself, but it does not quite work out...

def inverse_fourier(x, t):
"""Evaluate Fourier coefficients at point t.

Args:
    x: The Fourier coefficients.
    t: The time point at which to evaluate the function.
"""
x = x.flatten()
n = len(x)
k = np.arange(n)
y = x @ np.exp(2j * np.pi * t * k / n)

return y / n

plt.plot(np.array([inverse_fourier(fourier_coefficients, t) for t in range(50)]))

To give a bit more detail on the problem. I want to take the approximate difference between to signals that vary between different sinusoids over time, but are inherently phase shifted. So my idea was to take the last $k$ timepoints, estimate the Fourier coefficient, take the difference in Fourier space and evaluate the resulting function represented by the coefficients at time point $t$. Therefore, I want to know how I can evaluate the coefficients at $t$. I know the periods of the signals that I am varying between, so can go back enough $k$ points to capture all relevant frequencies.


texts: If $X$ and $Y$ are zero mean independent Gaussian random variables with different variances, what is the density of $\sqrt{X^2+Y^2}$

texts: Let $X\sim\mathcal{N}(0,\sigma_X^2)$ and $Y\sim\mathcal{N}(0,\sigma_Y^2)$ be independent Gaussian random variables. What will be PDF of $Z=\sqrt{X^2+Y^2}$ and $W=\arctan{\left(\frac{Y}{X}\right)}$. Will they still be Rayleigh and uniform distributions respectively?


texts: How to set the frequency cut-off for a low pass RC filter?

texts: How do I set the frequency for a low pass filter illustrated below?


texts: Circular Convolution as Cyclic Shift Operator

texts:
Given the following signal vectors:
$$ γ=[ψ_0,0,ψ_1,0,ψ_2,0,…,ψ_{N-1},0]^T\in \mathbb{R}^{2N}, ϕ=[1,\frac{1}{2},0,…,0,\frac{1}{2}]^T \in \mathbb{R}^{2N}$$
I want to show that the convolution of $γ$ and $ϕ$ is actually a cyclical-shift operator. i.e.:
$$γ*ϕ = γ + \frac{1}{2}T_1(γ) + \frac{1}{2}T_{2N-1}(γ) $$
Where $T_{t_0}$ is cyclical-shift operator with offset of $t_0$ places.
I tried to develop it according to the definition of (circular) convolution:
$$h_l=∑_{n=0}^{2N-1} γ_n ϕ_{l-n (mod \space 2N) } $$
However, I couldn't come up with something helpful...
How can I prove this identity? Thanks.


texts:
I know the FIR approach, I have seen IIR, to, but I'd like to know if it's possible to implement a Hilbert transformer in analog domain, i.e. with integrators instead of delays. Is it possible? If yes, how?

[edit]
After some searching on my own, I stumbled upon an answer in sci.electronics.design - "Analog delay??" (I'm not posting the link as it seems it involves some character string and I don't know what influence will that have). A simple search will do.

Still, this involves a translation of a digital delay into its analog counterpart, clock, etc, no continuous-time values. Could this mean that there is no other way to do it (analog)?


texts: Demodulating AM signal from carrier signal on Matlab

texts:
I'm trying to demodulate a voice signal (bandwidth approximately 5 kHz) from a carrier signal (100kHz).

  [y,Fs]=audioread('Signal1.wav');

  h=10*sinc(10*(-10:10));

  y1=filter(h,1,y);

But when I compare the signal with the "ground truth" signal provided, the size of the signal is 1/20th of the one I made. Did I miss some steps in the process? Thanks.


texts: Precise 5th and 7th harmonics of a sampled sine wave

texts:
Does anyone know in decibels (to 1/100th of a dB) what the theoretical 3rd, 5th and 7th harmonic of a 0dB fs 24-bit (i.e. full-level; 0dB = -8,388,607 to 8,388,607) sampled sine wave without dither will be? Where the harmonics are a product of the sampling process (caused by the effect of the discrete sampled steps). Better still, how I calculate the theoretical values for the harmonic components.

The fundamental is around 1,000Hz (happy to use exactly 1kHz). The sampling frequency is 48kHz. Asynchronous sampling.


texts: Find which signal was transferred in 32QAM

texts:
I've been assigned to solve a problem for 32QAM signals. The thing is, there are little to none sources in my native language to get any decent information on the topic.
The problem is (in its original form):

For 32QAM to determine which signal was transmitted using the minimum Euclidean distance between the symbols criteria, given that I and Q in the transmitted signal took values {+-1, +-3, +-5} and the received symbol sequence: [I, Q] = [1.5, -1], [-3.5, 4], [0.1, 0.1].

I assume, that [I, Q] can be replaced with [x, y] and the set of values is missing an A. If there are any misleading parts, I can specify them with my teacher, so you're welcome to ask.
It would be perfect if someone could explain the solution a bit for my further understanding or attach any sources to read.


texts: Cheaply generating coefficients for IIR/FIR given a cutoff frequency

texts:
I'm working on a DSP-based senior project, and we'd like to allow the user to provide a 20-20kHz cutoff specification for IIR/FIR high and low-pass filters on the fly.  This would be a simple problem on a desktop, but we're working on a TI 5535 DSP chip.  So we're looking for a C library for generating coefficients.  We're currently looking at adapting this code from C++ to C and cutting out the input stuff, but if there's a better solution we'd love to hear it.

For IIR - Any of Butterworth / Bessel / Chebyshev would work.  Whatever has the fastest runtime / lowest complexity is optimal.  Accuracy is not important (educational product, so as long as the output changes when they change the cutoff, we're fine).

For FIR - not really sure what to do here.  Raised cosine for the lowpass, ?? for the highpass.

We're sampling stereo audio at 44.1kHz/channel.  I'm unsure what other parameters we can set constant to make calculations simpler.  The number of poles can be held constant, and the rolloff and other specifications aren't terribly important to us.

The question - what are some existing methods or C libraries for generating filter coefficients from a cutoff frequency in the auditory range?


texts: Matlab: Difference between flat fading and Rayleigh multipath fading, terminologies and concept

texts:
The information that I put down are from the site http://www.gaussianwaves.com/2011/05/ebn0-vs-ber-for-bpsk-over-rayleigh-channel-and-awgn-channel-2/ . My QUestions are due to some misconception and will really appreciate if some help is provided to clear the concepts.

AWGN channel:

If $y$ is the received data and $n$ is the input, then signal is received  as $y = x + n$. This is an FIR model.

Rayleigh channel

The received signal here is $y = hx +n$ where $h$ can be $h_1$ or $h_2$. This is also an FIR model.

Case1: When $h = h_1$:

h_1 = 1/sqrt(2)*[randn(1,N) + j*randn(1,N)]; % Rayleigh channel creates the coefficients $h$ of length $N$ representing the channel, where $N$ is the number of data samples. This is a flat fading Rayleigh channel because the number of taps = 1. Please correct me if I am wrong in saying that the number of taps =1, so it is a flat fading. In general, number of taps = number of delays.

Case2: When $h=h_2$

If $h_2 = [1, 0, 0.5, -0.2]$ then also we call it Rayleigh but now it is also called multi-path. So, number of multipaths = length(h_2) = 4 ?

How is $h_2$ a Rayleigh random variable? What is the proper way to represent $h_2$ as rayleigh channel and the output $y$? Explanation with the help of code statement will be useful.
Theory says that there is Clarke's and Young's model for channel representation. Clarke's model uses FIR representation using sinusoids and cosines. Why do we need this model when we can simply apply the technique used to generate $h_1$ ?


texts: print a date and time in command window MATLAB

texts:
I use matlab R2015a, when I run my script I want print the date and the time of exection

I tried this code

Dtimes=datetime('now','TimeZone','local','Format','d-MMM-y HH:mm:ss Z');
% the first test
disp(' Date and time is ' Dtimes ); % But not work
% second test
disp(sprintf('date and time is %s  ',Dtimes)) %also not work

But with this command is worked :

disp(sprintf('date and time is %s  ',datestr(now)));

What is the problem?


texts:
In a book on optical solitons, the author says that a Gaussian signal pulse has the form-
$u(t) = \exp((\eta + i\beta)t^2)$
where $\beta$ is the chirp parameter.
My question is: what does the $i$ signify here? As a chirped signal is still real, what does the complex nature of the signal represent?


texts: Obtaining phase and magnitude of 2 spectral components after ADC

texts: I am currently implementing a digital filter to obtain the output of a PID controller. The control loop is carried out by applying a 1kHz sine wave at the input of the system and reading the output signal. The signal is distorted by the system and the output signal contains a 2nd harmonic component of the 1kHz input sine wave. I am currently applying FIR filters to obtain the phase and magnitude of the fundamental and harmonic tones but in order to reduce the noise I have to use very high order FIRs (around 1000 taps) and it is a very slow process. I would like to know if a matched filter would be a better approach, or if there are any other suggestions. Thanks!


texts: Ways to identify if music has vocal timbres in it?

texts:
The human voice is singular because it has a specific frequency spectrum and fast transitions between vowels and consonants.

What are good and computationally fast ways to detect if there is probably a voice contained in a noisy music, the amount of the song where there is vocals, and the type of vocals in different kinds of music, that perhaps have beats and guitars and other sounds mixed into the vocal?


texts: Can different Discrete-Time-Fourier-Series(DTFS) coefficients have the same discrete sequence in the time domain?

texts:
Please, check the following discrete periodic sequence when the period $N=2$.
$x[k]=\exp(j\frac{2\pi}{N}k),    N=\text{period}$
$..., x[0]= 1, x[1]= -1, x[2]= 1, x[3]= -1, ... , N=2$
According to my DTFS calculation, DTFS coefficient of above sequence is as follows.
$..., D_{0}=0, D_{1}=1, D_{2}=0, D_{3}=1, ... , N=2 $
If I shift above DTFS by 1 in the frequency domain like this,
$..., D_{0}=1, D_{1}=0, D_{2}=1, D_{3}=0, ...$
Inverse DTFS of above coefficient is as follows according to my calculation.
$..., x[0]=1, x[1]=1, x[2]=1, x[3]=1, ...    N=2$
Actually, if DTFS of $x[k]$ is $D_{n}$, then DFTS of $\exp(ja\frac{2\pi}{N}k)x[k]$
will be $D_{n-a}$ by DTFS frequency shifting definition.
So, DFTS of $\exp(j(-1)\frac{2\pi}{N}k)x[k]$ will be $D_{n+1}$.
In this case, $\exp(j(-1)\frac{2\pi}{2}k)x[k]$ is
$\exp(-j\frac{2\pi}{2}k)$ $\exp(j\frac{2\pi}{2}k)=1$.
Shifting by $-1$ shows the same result in this case because $N=2$.
So, I think above result is correct.
I know DTFS of $x[k]=1$ with $N=1$ is as follows.
$..., x[0]=1, x[1]=1, x[2]=1, x[3]=1, ...    N=1$
My conclusion is as follows.
$x[k]=..., 1, -1,  1, -1,...    \longleftrightarrow  D_{n}=..., 0, 1, 0, 1,...   N=2$  case 1
$x[k]=..., 1,  1,  1,  1,...    \longleftrightarrow  D_{n}=..., 1, 0, 1, 0,...   N=2$  case 2
$x[k]=..., 1,  1,  1,  1,...    \longleftrightarrow  D_{n}=..., 1, 1, 1, 1,...   N=1$  case 3
I can't understand the last 2 lines.
I think case 2 and case 3 should have
different sequences in the time domian.
What did I miss?


texts: Are the DFT/FFT multiplications in the fourier domain complex multiplications?

texts:
If you want to compute a convolution with a DFT/FFT, you will have to do a DFT on the image and a DFT on the kernel.

And you will need to multiply the spectrum of the image with the spectrum of the kernel.

Is it correct that the  multiplications in the fourier domain are complex multiplications? So, that those point-wise multiplications will require 4 real multiplications and 2 real additions?


texts: Locking on to a square wave signal with minimum oversampling

texts:
I'm designing a device that will have an IR photodiode connected to a low power microcontroller's ADC pin. At times, another device will be transmitting a 48KHz square wave, and I'd like to be able to detect when that is the case entirely in software, as well as being able to detect individual high and low transitions.

My understanding is that a PLL is the correct solution to this problem. However, I have extremely limited DSP knowledge, so I have a fuzzy idea at best on how to actually implement this in software, and much of the documentation available goes well over my head with a lot of domain specific terminology. This is further complicated by the low power nature of the microcontroller - I have limited processor cycles available, and thus I'd like to sample as infrequently as possible while still reliably locking on to the signal. If it helps, the sampling timer can be adjusted to match the observed frequency of the input signal.

Can anyone suggest an approach to this that is as straightforward and understandable as possible, while also minimizing CPU resources consumed?


texts: How to add noise to time series data in Python?

texts: I have time series data with $x$, $y$, $z$ sampled at $50\textrm{ Hz}$ from a sensor and I would like to add a uniform and Laplace noise to it. How would I achieve this in Python?


texts: Is there processing gain for FMCW using heterodyne-style receiver as opposed to matched filter?

texts:
Beat signal of a single target will be a sinusoid in the idealized world, so theoretically the signal processing gain of an FMCW pulse correlated with Tx waveform in this way should be analogous to FFT processing gain of a single tone, right?
Just trying to wrap my head around what to expect in the field when playing with waveform parameters and what kind of SNR before/after I can expect to see. Is it really just a product of the duration of wave form and bandwidth? Looking into buying a textbook on the subject of FMCW so if you have any recommendations feel free to drop em here


texts: Quantifying the energy of a voiced frame

texts:
What should the energy, pitch frequency, and zero crossing rate of a voice/unvoiced/silence frame in signal be?

My professor said that it is experimental, but I do not have any idea what the factors of these are and where I should begin from.

I would appreciate any guesses based on experimental data or any limits on the quantities.

I've tried a zero crossing rate between 2000 and 94000, energy between 0.0001 and 0.7, and
pitch frequency between 180 and 500.


texts: How to detect that the same word has been spoken?

texts:
Say, I pronounce a word and record it using a certain microphone. This record is treated as a master word. Then I record the same word again using the same microphone. I need to somehow compare the second sample to the master word and tell if it's the same word or not. I was thinking about a simple comparison with a threshold in percent. Say, > 80 % match means that it's most likely the same word, < 80 % match suggests that it's highly likely a different word.
It's important that this in no way should be any kind of speech recognition system, that would definitely be too complicated and an overkill in my situation. I'd rather call it a sound matching system
How should I go about this and what programming language/framework/library can work best for my goal?


texts: Difference between $\tt freqz()$ and $\tt plot(abs(fft()))$ in MATLAB

texts:
I've plotted the response of a high-pass filter x=[1 -1];

freqz(x) gave me

where as the magnitude plot for the fft gave me

Why are the two plots different? Both are frequency responses, then why are my plots are diffferent. What does freqz() do?


texts: Phantom harmonics when using cosine windows why do they appear and how to avoid them?

texts:
Given an $L$ order cosine window, it is possible to show that the width of the main lobe is given by:
$$\omega_w = \frac{2 \pi L}{(2N+1)}$$
Where $L$ is the order of the window, $N$ is the maximum harmonic in the signal and therefore $2N+1$ is the number of samples taken over a period. If we consider the relationship
$$\Omega = w T_s$$
then the width of the window in the frequency domain is given by
$$\Omega_w = \frac{2 \pi L}{(2N+1)T_s}=\frac{2 \pi L}{T_0}$$
Where $Ts$ is the sampling time and $T_0$ is the observation interval.
The last formula clearly shows that it is possible to reduce the width of the main lobe either by oversampling by an oversampling ratio of $L$ or by increasing the observation interval by a factor of $L$.
Let's assume the following:

Sampled signal is a pure sine wave with frequency $f = \frac{1}{T}$. Then $N=1$.
The signal is sampled according to the synchronous sampling condition ($(2N+1)T_s = T_0 = T$) so no leakage error occurs.
The Hanning window is used. The width of the main lobe should be $\Omega_w = \frac{4 \pi}{(2N+1)T_s}$ since $L = 2$.

Now, if we call $P(j\Omega)$ the spectrum of the sampled signal, $P(j\Omega) = W(j\Omega) \star S(j\Omega)$ where $W$ and $S$ are the Fourier transforms of the window and of the signal, respectively and $\star$ is the convolution operator.
What would $P(j\Omega)$ look like? In the book I use, it looks like this:

STATEMENT A:
In my understanding at least, only $c$ is real (only $c$ is a replica of the real spectrum of the signal) while $a$ and $b$ are a byproduct of the convolution with the window due to its main lobe's width. In general, the book says, it is possible to say that if we are using an $L$ order cosine window, $2L-2$ phantom (or fictitious if you will) components will appear.

STATEMENT B:
If we increase $T_0$ by $L$ then the width of the Hanning window becomes the same as the one of the rectangular window ($\Omega_w = \frac{2 \pi}{(2N+1)T_s}$), therefore the fictitious harmonic components should not appear. But isn't this the same as increasing the number of samples by a ratio $L$?

I would like to understand if the two above statements are correct, in particular if the origin of the fictitious harmonic components is indeed the convolution and how one can avoid these components.
The book that I mention is Digital Signal Processing for Measurement Systems by D'Antona and Ferrero.


texts: Multiple choice on sampling and aliasing

texts:
I found some multiple choices in a well known book . The problem is that I don't get the answers in some and I want to do so.
Question 1:
The signal x(t) with Fourier transform $X(j\omega) = u(\omega)- u(\omega -\omega_0 )$ can undergo
impulse-train sampling without aliasing, provided that the sampling period $T <\frac{2\pi}{\omega_0}$.
Answer: True
Okay here's what I think: $x(t)$ should be sampled with at least:
\begin{align}
ωs &> 2B \\
&> 2(B_1-B_2)\\
&>2(\omega_0 - 0) \\
&> 2\omega_0 \\
\frac{2\pi}{T_s} &> 2\omega_0 \\
 T_s&<\frac{\pi}{\omega_0}
\end{align}
What do I miss?
Question 2:

Answer: $50\pi$
For $x(t)$, $\omega_{\text{max}} = 50π + 50π = 100π$ and so $ω_s > 2\omega_{\text{max}} = 200π$.
So , if we actually sample with $\omega_s<\omega_{\text{Nyquist}} = 200π$ , then there is no chance to avoid aliasing for all $X(j\omega)$ spectrum but we can gain a proper copy of it as long as $\omega_0 < \frac{\omega_s}{2} = 75π$. Again what do I miss?


texts: This result has been used in the second last line of the pic. I don't know why it's true. Both functions are zero for $t$ not equal to $-b/a$. But at $t=-b/a$, a scaling factor $1/|a|$ has been introduced in the second function.


texts:
I am trying to implement a Weighted Matching Pursuit algorithm. Essentially I would like to weight the atoms of a Dictionary (redundant matrix) with a frequency masking function, like the one specified for MP3. I have read some of the papers, and theory of frequency masking, but I am still having trouble figuring out how to generate the function in Matlab. My dictionary contains frequencies ranging from -fs/2 to fs/2. fs = sampling frequency.

If anyone can help me out with a quick way to generate this kind of function I would be most grateful.


texts: Similarity in signal plots but dissimilarity in sensed content

texts:
I have two audio signals that look incredibly similar when plotted. However, when I play them, they sound quite different. The only difference I can even see in their plots is that one signal looks darker than the other. Why would this be happening?

FYI: The signal on the top is one that has been reconstructed from a reconstruction algorithm and the one in the bottom is the original one.

Here's some further data about how the signal has been produced. The original signal is a sample audio data I downloaded off the internet. The purpose of my experiment is to produce a close enough signal from only the amplitude spectrum and without the phase information. For similar, but slower algorithms, one can look at the Griffin and Lim algorithm and also the single pass spectrogram inversion paper.

In my case, I'm getting a waveform very close to the original signal, but there's a significant "noise" when listening to the output waveform.

So in short, the question is: waveforms looking so close to each other sound very different which is not very intuitive. Also, when signals look so similar, it is difficult to discern that there's any kind of noise involved. So what kind of hidden noise is this?

One last piece of information: the correlation function gives 0.69 as the output when comparing these two signals. Not sure if that's useful


texts:
I have created a little classifier that listens for certain keywords. Sometimes it will not react to a keyword, but if I shout the word, then the classifier will be able to pick it up.
The classifier uses MFCC as features, which is basically a Short-time Fourier Transform with some extra stuff. This tells me that when I shout, the frequencies that the classifier are detecting are actually different from the non-shout.
So I'm wondering: does the volume of my shout change the frequencies in my speech, or am I somehow enunciating differently when I shout?


texts: Optical Flow: Method to estimate a very small displacement

texts:
As I was searching for available resources( here, here, here and many more), I realize that many developed optical algorithm is to estimates LARGE displacement. However, there is very limited information or journals that put effort on estimate very small displacement. The small displacement I'm taking here is at sub-pixel level. The smaller the resolution the better.

So my question, you guys have any idea in mind that might inspire a new method, or suggest an existed method (that I might not be aware of/that can be improved)? Any method/suggestion/comment is very welcomed!


texts:
I know that in many edge detection algorithm the first step is choosing a scale to work on.
We actually operate on the chosen scale by applying a gaussian blur on the original image.
I wanted to know , what is the relation between the detected edge's width to the sigma.

e.g. If I want to detect lines 5 pixel wide , which sigma should I use?
How should I change it if I want to find a line 10 pixels wide?

Thanks in advance,


texts: Compressing a time-domain signal in the frequency domain Matlab

texts:
I'm new to frequency domain analysis. Let's say I have a frequency domain representation from a time domain signal. What I want to do is taking frequencies with a low amplitude out of this frequency domain representation and then convert this frequency domain representation back to the time domain such that the result is almost the original time domain representation.

In Matlab what I want to do is something like this:

timeDomain(1:50) -> apply fft(timeDomain) -> inspect frequency domain and filter out low amplitude frequencies yielding a frequency domain vector x with length(x) probably (way) smaller than 50. (so x should contain only the frequencies with the highest amplitudes of the signal) -> apply ifft on this filtered complex vector -> Get an approximation of the original time domain signal (timeDomain(1:50) that has also length 50.

So my aim is to represent the time domain signal timeDomain(1:50) with less dimensions in the frequency domain without too much information about the signal getting lost, thus being able to reconstruct the time domain signal, so basically compression.

Thanks for your time.


texts:
I am trying to calculate shannon entropy of CWT. I am not sure if I am doing it right.
Assume that $W(a_i,t), i=1;2;...;M$ is a set of wavelet coefficients. The Shannon wavelet entropy is calculated by:
$E=-\sum_{i=1}^{M}d_i log(d_i)$ $\rightarrow$   where $d_i=\frac{|W(a_i,t)|}{\sum_{j=1}^{M}W(a_j,t)}$
I am confused how to calculate $E$. for example I have a coefficient matrix with size of $M\times N$, $M$ is scales number and $N$  is time segments. first I have to calulate $d_i$, this is my main problem. this is the wavelet coefficient matrix :
 $W_{M\times N} =
 \begin{pmatrix}
  w_{a_1,1} & w_{a_1,2} & \cdots & w_{a_1,N} \\
  w_{a_2,1} & w_{a_2,2} & \cdots & w_{a_2,N} \\
  \vdots  & \vdots  & \ddots & \vdots  \\
  w_{a_M,1} & w_{a_M,2} & \cdots & w_{a_M,N}
 \end{pmatrix}$

hmm i am pretty sure i am wrong, can anyone help?
for example tell me how can i calculate $d_4$?
Here I have right a little Matlab script to calculate shannon entropy of CWT.
Is it right or wrong? and what should I do?

 [M,N]=size(coeffs);
for js=1:M
Ej(js)=sum(abs(coeffs(js,:)));
end;
Etot=sum(Ej);
Pj=Ej./Etot;
%shannon entropy
shan_entr=-sum(Pj.*log(Pj));


texts: How to Determine Covariance Matrix $Q$ and $R$ in Kalman Filter

texts:
I am implementing getting orientation from smartphone. I want to use Kalman filter and should determine process noise covariance matrix $Q$ and measurement noise covariance matrix $R$. (newbie to Kalman filter)

I don't have any idea how to determine $Q$. What I think about $R$ is as follows:

state vector : quaternion from (accelerometer + gyroscope)

(1) My phone is stand still. I get covaraince matrix from Matlab

1.0e-04 *

0.0000    0.0005    0.0035   -0.0000
0.0005    0.0063    0.0411   -0.0002
0.0035    0.0411    0.2881   -0.0014

-0.0000   -0.0002   -0.0014    0.0000

(2) My phone had been moved for 5 seconds.

covariance matrix is

0.0417   -0.0533   -0.0008   -0.0014
-0.0533    0.0784    0.0015    0.0018
-0.0008    0.0015    0.0001    0.0001
-0.0014    0.0018    0.0001    0.0001

Is there anyone to help?

(added)

Details are omitted.

case 1: Kalman Filter

The row data from my phone is $p, q, r$ (angular velocity). I omit the conversion equation between angular velocity and quaternion.
\begin{align*}
x_{k+1} &= Ax_k+w_k \\
z_k &= Hx_k + \nu_k \\
Q &: \text{ covariance matrix for }w_k\\
R &: \text{ covariance matrix for }\nu_k
\end{align*}

$$ \begin{bmatrix}\dot q_1\\\dot q_2\\\dot q_3\\\dot q_4\end{bmatrix}=\frac 12 \begin{bmatrix} 0&-p&-q&-r\\p&0&r&-q\\q&-r&0&p\\r&q&-p&0 \end{bmatrix}
\begin{bmatrix} q_1\\ q_2\\ q_3\\ q_4\end{bmatrix}$$

$$ \underbrace{\begin{bmatrix} q_1\\ q_2\\ q_3\\ q_4\end{bmatrix}_{k+1}}_{x_{k+1}}=\underbrace{\left( I + \Delta t\cdot \frac 12 \begin{bmatrix} 0&-p&-q&-r\\p&0&r&-q\\q&-r&0&p\\r&q&-p&0 \end{bmatrix}\right)}_{A}
\underbrace{\begin{bmatrix} q_1\\ q_2\\ q_3\\ q_4\end{bmatrix}_k}_{x_k} $$
$$ H=I$$
My guess for covariance matrix is as follows: (but I don' know how to infer..)
$$ Q = 0.001I, \quad R=10I.$$

case 2: Extended Kalman Filter

\begin{align*}x_{k+1} &= f(x_k) + w_k \\
z_k &= h(x_k) + \nu_k\\
Q &: \text{ covariance matrix for }w_k\\
R &: \text{ covariance matrix for }\nu_k
\end{align*}
$$ A = \left.\frac{\partial f}{\partial x}\right|_{x_k} ,\quad H = \left.\frac{\partial h}{\partial x}\right|_{x_k} $$
\begin{align*} \begin{bmatrix} \dot \phi\\ \dot \theta\\ \dot \varphi \end{bmatrix}&=
\begin{bmatrix} 1&\sin\phi\tan\theta & \cos\phi\tan\theta \\
0&\cos\phi & -\sin\phi \\
0&\sin\phi\sec\theta & \cos\phi\sec\theta \end{bmatrix} \begin{bmatrix}p\\ q\\  r \end{bmatrix} \\
&= \begin{bmatrix} p+q\sin\phi\tan\theta+r\cos\phi\tan\theta \\
q\cos\phi-r\sin\phi \\
q\sin\phi\sec\theta + r\cos\phi\sec\theta \end{bmatrix} \\
&= f(x) + w
\end{align*}
$$z = \begin{bmatrix}1&0&0\\0&1&0\end{bmatrix}\begin{bmatrix}\phi\\\theta\\\varphi\end{bmatrix}+\nu = Hx + \nu $$
$$A = \begin{bmatrix} \frac{\partial f_1}{\partial\phi} & \frac{\partial f_1}{\partial\theta} & \frac{\partial f_1}{\partial\varphi} \\
\frac{\partial f_2}{\partial\phi} & \frac{\partial f_2}{\partial\theta} & \frac{\partial f_2}{\partial\varphi} \\
\frac{\partial f_3}{\partial\phi} & \frac{\partial f_3}{\partial\theta} & \frac{\partial f_3}{\partial\varphi}
\end{bmatrix} $$

(I emphasize that details are omitted.) In this case, also I don't know how to infere $Q, R$.


texts: Power Spectral Density computation and units

texts:
I want to make some calculs of power spectral densité of signal. For example a real voltage signal (physical unit : $V$) in time $g(t)$, its fourier transform $G(f)$ and $S_g(f)$.
   As far as I know, the power spectral density units for g(t) is $V^2/Hz$. However, I find in various sources the following equation to calculate the power spectral density :
  $S_g(f)=|G(f)|^2=G(f)G(-f)$

I do not understand this "definition" (of course the most fundamental definition of PSD is given by the Wienner-Kintchine theorem) , as it leads to a spectral density in $V^2/Hz^2$.

Regards,
      Mike


texts: How to Use Maximum a Posteriori Probability (MAP) in Classification Task

texts:
I have a 2D image defined on a region $\Omega$. Let $I: \Omega \to R$ be a gray image. Assume that the region can be separated into $N$ sub-regions $\Omega_i$ such that $$\forall i,j=1,\ldots ,N:\Omega_i \cap\Omega_j=\emptyset$$ and $$\bigcup_{n=1}^{N}\Omega_n =\Omega$$
To simplify, we assume $N=3$ and it shows in below figure.

We denotes $p(x \in \Omega_i|I(x))$ is posteriori probability of region $\Omega_i$ given $I(x)$. Based on the Bayesian rules, we have
$$p\left(x \in \Omega_i|I(x)\right)=\frac{p\left(I(x)|x\in \Omega_i\right)p(x\in \Omega_i)}{p(I(x))}$$
The segmentation try to assign pixel $x$ to each region based on posteriori probability. That goal is that finds maximum a posteriori probability (MAP) of pixel $x$ in each region. It can express by formula
$$\left\{ \Omega_{i,x} \right\}_{i=1}^{N}=\arg \max_{\Omega_i} \left\{\prod_{i=1}^{N} p(I(x)|x \in \Omega_i)p(x \in \Omega_i )\right\} \tag{*}$$
Assuming that the pixels within each region are independent , the MAP will be achieved over all image domain $\Omega$
$$\left\{ \Omega_i \right\}_{i=1}^{N}=\arg \max_{\Omega_i} \left\{\prod_{i=1}^{N}\prod_{x \in \Omega_i} p(I(x)|x \in \Omega_i)p(x \in \Omega_i )\right\} \tag{**}$$

Are my formulas $(*)$ and $(**)$ correct?
If not, Could you help me to correct them?

The reference is Geodesic Active Regions and Level Set Methods for Supervised Texture Segmentation.


texts: When is cubic spline interpolation better than an interpolating polynomial?

texts:
The following plot is a slight variation of an example in a text book. The author used this example to illustrate that an interpolating polynomial over equally spaced samples has large oscillations near the ends of the interpolating interval. Of course cubic spline interpolation gives a good approximation over the whole interval. For years, I thought high order polynomial interpolation over equally spaced samples should be avoided for the reason illustrated here.

However, I recently found many examples of bandlimited signals where a high order interpolating polynomial gives less approximation error than cubic-spline interpolation. Typically an Interpolating polynomial is more accurate over the entire interpolating interval when the sample rate is sufficiently high. This seems to hold when the samples are equally spaced with a sample rate at least 3 times greater than the Nyquist frequency of the signal. Furthermore, the advantage over cubic spline interpolation improves as (sample rate)/(Nyquist frequency) increasees.

As an example, I compare cubic-spline interpolation with an interpolating polynomial for a sine wave with a Nyquist frequency of 2 Hz, and a sample rate of 6.5 Hz. Between the sample points, tthe interpolating polynomial looks exactly the same as the actual signal.

Below I compare the error in the two approximations. As with the first example, the polynomial interpolation does worst near the beginning and end of the sample interval.  However, the interpolating polynomial has less error than a cubic spline over the whole sample interval. The interpolating polynomial also has less error when extrapolating over a small interval. Did I discover a well known fact? If so, where can I read about it?


texts: Suppose we have a system whose impulse response h has length K and fed with an input x that has length N. Then it is known that the output y has length M = K + N -1. This shows us the convolution matrix which relates x and y has size M x N. By existence theorem, the problems in the format $A_{m x n}x = b$ has a solution if and only if $m \leq n$. Since $m > n$ always holds, it seems to me that there will be always no solution for FIR filter deconvolution problem. Is this correct becase in a somewhere, I read about the possibility of existence.


texts: Is it possible to use k-means algorithm with just one vector

texts:
Suppose I have a vector $X = (x_1, x_2 , . . . ,x_n)$, $x_i$ is the maximum of $X$ and $x_k$ is the minimum. Is it possible to use k-means algorithm to cluster the values in vector $X$ into two cluster. First, the values which are near to maximum value and the other cluster includes the values which are near to minimum value?
If not, Is there another algorithm can do that task?

Thank  you


texts: Energy Detection in Presence of Colored Gaussian Noise

texts:
Before asking my question, let me introduce the context:

For spectrum sensing based on energy detection, which has been widely studied in presence of AWGN, the optimal detection threshold is computed as:

$$ \lambda = 2~\sigma^2_\omega ~ \sqrt{N}~ (Q^{-1}(P_{fa}) + \sqrt{N})~~~~[1]$$

where $\sigma^2_\omega$ is the noise variance, and $P_{fa}$ is the design false-alarm rate (false-alarm probability).

The detection algorithm works as follows:

We take an interval of $N$ samples of the received signal $x[n] = k~s[n] + \omega[n]$, where $k = 1$ if the signal of interest is present and $k = 0$ otherwise. Here both $s[n]$ and $\omega[n]$ come from zero-mean Gaussian wide-sense stationary processes.

The received energy is then $\Lambda = \sum_{n = 0}^{N-1} |x[n]|^2$, and $\Lambda$ follows a $\chi^2$ distribution.

Now, regarding the presence or not of a signal of interest during the considered interval, it is assumed that $k = 1$ if $\Lambda > \lambda$, and $k = 0$ otherwise.

Of course, equation [1] is not optimal in practical scenarios, since AWGN is just a model not possible in real life. However, let's just accept it as a good approximation for wide-band noise (nearly-uncorrelated noise). After applying an LTI pass-band filter to such a noise, I would obtain Gaussian colored noise (no longer wide-band), and thus uncorrelated samples is no longer a good approximation. As a consequence, not much can be said about the distribution of $\Lambda$, due to the introduced correlation in $x[n]$ after filtering.

My question is:
Is it still justified to assume that the detection threshold $\lambda$ is linearly proportional to the noise variance (after the pass-band filter)?

My results in this regard: A practical application on acoustic signals I am working on suggests a non-linear relation between $\lambda$ and the noise variance, once fixed a false-alarm probability.

I mean, I fix a threshold, then I run the application in absence of a signal of interest, and after a long run I measure the false-alarm rate. By increasing the noise power, for example 10 dB, and then doing the same for the threshold, I obtain a false-alarm rate quite different from the previous run. However, I don't see a proper justification for this behavior other than a non-linear relation between $\lambda$ and the noise variance. And in such a case, I am unable to explain why is this so.

I am using as the receiver a mobile (its microphone), and as the transmitter a speaker.

Should I empirically determine what is this non-linear relation, or is there any theoretical result in this regard that I could use?


texts: Sampling rate conversion and polyphase decomposition

texts:
Can someone please explain the answer to this problem:

The sampling rate of a baseband signal shall be converted from 5 kHz to 7.5 kHz using the efficient polyphase implementation with a 'rotating' switch. How many polyphase components are required and how many steps must the switch be rotated?


texts: Mean delay for a M tap harmonically weighted filter

texts:
I am trying to derive the mean delay for $M$ filter taps harmonically weighted:

$$y[n] = \frac{M\,x[n] + (M-1)x[n−1]... + 1 \cdot x[n−M+1]}{\tfrac12 M (M + 1)}$$.

For a uniformly weighted filter, I used group delay formula and solved it but I am not able to solve this question in a similar way. If I find frequency response $H(e^{j\omega})$ for the above mentioned equation, I end with a arithmetic-geometric progression

$$H(e^{j\omega})=M+(M-1)e^{-j\omega}+...e^{-j\omega(M-1)}$$.

I don't know whether my approach is right or not.

Also, if this is a non-linear phase filter, how should I solve it?

Thanks in advance.


texts:
I came across various works in the literature of document analysis addressing the problem of detecting base lines in multi-oriented texts. But no amount of search could give me satisfactory answer to the following quest.

Is there an algorithm to detect the base line which neither use a histogram based method or any kind of threshold parameter?


texts: Determining the Filter Coefficients of an FIR Filter

texts:
I need to find the filter coefficients of an FIR filter that will block sinusoids of frequency $200\ \rm Hz$ if the sinusoid is sampled at $1.2\ \rm kHz$.

I feel like this is a fairly simple problem, but I'm not exactly sure how to go about it.

I am thinking I need to convert the frequencies to rad/s, then use these values in the transfer function, but am not sure how to find this transfer function.

EDIT:

I have to choose from one of the possible options:

$$ \{1, 1, 1\} $$
$$ \{1, -1, 1\} $$
$$ \{1, 0, 1\} $$
$$ \{1, \sqrt{2}, 1\} $$

or None of the Above.

From the coefficients above, I know the following equations can be created:

\begin{align}
H\left(e^{j\omega}\right)&=1+e^{-j\omega}+e^{-2j\omega}\\
H\left(e^{j\omega}\right)&=1-e^{-j\omega}+e^{-2j\omega}\\
H\left(e^{j\omega}\right)&=1+e^{-2j\omega}\\
H\left(e^{j\omega}\right)&=1+\sqrt{2}e^{-j\omega}+e^{-2j\omega}
\end{align}


texts:
I'm receiving PCM samples through stream which is one channel, 8000 samples per second. I figured out way of finding out loudness by calculating 'force' of signal by adding abs value of every next 1000 samples and then dividing by samples count which is 1000. It works somewhat good but sometimes my calculations are going crazy, my sum of absolute values of samples seems to constantly build up and reach ridiculous values... so I think that I'm doing something wrong there. Funny thing is that when calculations are going crazy everything will go back to normal when I hit 'sound input device' with my hand ;)

Anyway. My question is: how to properly calculate loudness of set of PCM samples?


texts: How to develop a spectrogram (2D array) from audio signal?

texts:
I have developed a spectrogram in Python using Scipy.Signal.Spectrogram. But I need a complete understanding of data. Here I am not asking about plotting and color selection etc. I am more into data (numbers). I am attaching a picture, please have a look :

In image 1, you can see that I have an array of frequency mapped to [0, Fs/2] i.e [0,24000] in my case.

In image 2, Time is mapped from [0,10sec] and the total array length is 2141.

In image 3, Spectrogram has been computed in a 2D array.

I want a clear understanding that how these array of Frequencies and Times have been developed. What is the data that the spectrogram holds in the 2D array? Is it a log magnitude of frequency domain components?

I need some clear steps as I want data in 2D array as it can be seen in the 3rd Image. Here I am not into plotting thing but more into a clear understanding of data behind the spectrogram. I want to get this same data in C++.