texts: Formula to calculate Cepstral coefficients (not MFCC)

texts:
I am currently working with non-audio signals of which I would like to calculate the Cepstrum coefficients with Python so that I can use them with machine learning algorithms. That's probably a quite basic question, but somehow I'm struggling at this point.

I already found some packages in Python that can be used to calculate the MFCCs. However, since these refer to the Mel Frequency, they are unfortunately not suitable for my application.

I already found out that the cepstrum of the signal can be calculated in Python as follows (from this website):

powerspectrum = np.abs(np.fft.fft(signal))**2
cepstrum = np.fft.ifft(np.log(powerspectrum))

That's working so far. But with which formula do I get the (n first) cepstral coefficients?

(I know that I should apply a window function on the signal first - my question is mainly about the approach to calculate the cepstral coefficients)


texts: Why is my signal attenuated after applying a HighPass-filter?

texts:
I am new to DSP, but I've been attempting apply a butterworth high-pass filter on a 1D array (Values peeking at 300~) to eliminate any noise, yet when applying a high pass filter, I'm ending up with an attenuated result?

So far I have attempted to approach the problem with different specs, but it seems that my lack of experience in Signal Processing is holding me back.

def butter_highpass(highcut, fs, order):
    nyq = 0.5 * fs
    high = highcut / nyq
    b, a = butter(order, high, btype="highpass")
    return b, a

fs      = 30.0            #Sampling rate Hz
b, a = butter_highpass(5, fs, 5)
filtered      = signal.lfilter(b, a, data)

Result:

My Goal : Reduce all noise to 0 while retaining the peak values.

Any suggestions or recommendations are more than welcome


texts: Second-order sections from complex coefficients

texts:
I have an impulse response ir that I need to approximate with an iir because of realtime requirements. For stability I want to break down the 10th order IIR generated with Steiglitz-Mcbride to biquads / second-order sections. However, I can't use tf2sos with the complex transfer function that stmcb() generates like the following:

[b,a] = stmcb(ir,10,10);
[sos,g] = tf2sos(real(b),real(a));

When I compare bode plots of

freqz(b,a)

and

freqz(sos)

I find that the contour is preserved but there is a large difference in magnitude and phase. Am I missing some scaling factor or limitation? Thanks for any help.


texts: What is the purpose of the recurrence relation in low and high pass (audio) filters?

texts:
Newbie here. I'm familiar with time and frequency domains from math. I think I understand what low and high pass filters do, sonically. But digging into an implementation of tone in csounds, I think I'm missing something fundamental. Take this description:

A first-order recursive low-pass filter with variable frequency response.
tone is a 1 term IIR filter. Its formula is:
yn = c1 * xn + c2 * yn-1
where

b = 2 - cos(2 π hp/sr);

c2 = b - sqrt(b2 - 1.0)

c1 = 1 - c2

I see that there is a recurrence relation. But I'm reading about IIRs and filters in general I don't see any explanation of why it's there or what it does.

Would performing a single pass of this formula not be considered performing a "true" low or high pass filter?

Or is that it takes multiple passes (infinite, in theory) to completely "remove" the frequencies desiring to be removed, so it's a tradeoff between iterations versus computation?

I'm gravitating toward the latter, but I don't see where that is being said. Can someone point me to the essential statement of it on Wikipedia or somewhere?


texts:
lets say I have 3 audio filters for a realtime DSP application. For simplictity each has length 256 (as well as the input signal). The filters should work in series. Starting with filter IR h1(n) and input x(n), zeropadding to 512, FFT, then multiplication and IFFT, delivers the first stage output Y1(n).

Since they input vector was padded to 512 the output has 512, too.

Second filter stage: input to this stage is y1(n) which has already 512 samples, and has to be convolved with filter h2(L=256). So y1 and h2 have to be zeropadded to 788, FFT, Multiplication, IFFT... Output y2 has lenght 768

Third stage: will add another padding, so final output vector will have length 1024

Since my app is realtime I would, for performance reason, like to avoid multiple FFT->iFFT for each stage. Moreover there will be a growing input length for the FFT with each filter stage.
I would prefer doing FFT for each filter individually to get H1,H2,H3 then do X*H1*H2*H3, then IFFT. As far as I understood from a previous post, this would work, if I would zero pad each filter (and input) already in the beginning to length L=N+N1+N2+N3. Unfortunately this I want to avoid since some filters are laying already pre-FFt´ed in big buffers already from the app start to grant fast access, and others (Equalizers e.g.) can be optionally switched on/off during runtime. So I have no way to know at program start to know how much "padded" zeros I will need to apply to those pre-FFT´ed filters.

Is there a way to get arround this? Maybe zero padding those static filters in frequency domain "on the fly" if there should be added e.g. another EQ filter on runtime? If yes, how?


texts: How to quickly calculate all the line segments when their coordinate in the Hough transform is given

texts:
Hough transform is a very powerful tool to find dominant straight lines in a black-and-white image. A very common practice of analyzing Hough transform result is to find some local maximum points in the Hough transform image, and then try to analyze the straight line segments in the image domain that correspond to the local maximum point. Then the question arises: how can we quickly calculate all the line segments when their coordinate in the Hough transform is given? For the time being, I have a awkward way to solve this problem:

first calculate the straight line function that correspond to the
Hough transform point.
after that, find all the line segments    that have intersections
with the calculated straight line. In the    meantime, I make sure
that all the pixels in the line segments are    near to the
calculated straight line by setting a some threshold.

Any ideas will be appreciated.


texts: I saw a lot of use of tensor in signal processing. What is the intuition behind it? Is it simply a common representation for audio (1D), image (2D) and video (3D) signals?


texts: Coherent averaging of polyphase components of a signal

texts:
I want to subsample a multi-band signal to introduce controlled aliasing. My frequency components are uniformly distributed across the spectrum (comb structure) and occupy a relatively narrow bandwidth each. The spectrum can be said to be sparse.

I achieve decimation through picking every $n$-th sample from the original signal. This of course does not take advantage of the remaining samples, which are simply rejected. By adding a phase (sample) shift in my decimation routine, I can get $n$ under-sampled signals in total, which are known as polyphase components. I can also incoherently average them (periodogram averaging) in the spectral domain, since the spectrum is of my interest rather then the time-domain representation.

Is it possible to coherently (in phase) average such undersampled polyphase components of my signal? What phase shift should be applied after down-sampling with a factor $n$? In other words, can we simultaneously undersample (decimate) and oversample (average) to take advantage of the remaining samples for quantization noise suppression?

A typical scenario would be sampling a 1 GHz wide comb signal with a spacing of 10 MHz between the teeth, which we can decimate 10 times to have a spacing of 1 MHz and a bandwidth of 100 MHz. If we sample the signal with an 8-bit analog-to-digital converter, can we use the remaining 9 polyphase components to enhance the vertical resolution of the signal after subsampling?

Thank you,
Lukasz S.


texts: How to represent impulse function in 2D?

texts:
To be more specific I want to show that impulse function in 2D can be represented as $β(r)=δ(r)/πr$.

Also I want to show that each projection of a two dimensional impulse function at the origin is a Delta function.

I could not find any useful information about these two problems, if anyone could help I will be appreciated.


texts:
I have a signal with multiple frequencies, and two of them, one of which is my main frequency, overlap.

Are there any techniques that could separate two frequencies that almost overlap?
I can increase my sampling rate but not my window size

I have already tried SVD,PCA,ICA but none of the methods managed to decompose those two frequencies.

N=200
x = np.zeros(N)
y = np.zeros(N)
z = np.zeros(N)
k = np.zeros(N)
x[0] = np.cos(2*np.pi)*np.cos(2*np.pi*0.1)
y[0] = np.cos(2*np.pi)*np.sin(2*np.pi*0.1)
z[0] = np.sin(2*np.pi)*np.cos(2*np.pi*0.1)
k[0] = np.sin(2*np.pi)*np.sin(2*np.pi*0.1)
for i in range(1,N):
  x0=x[i-1]
  y0=y[i-1]
  z0=z[i-1]
  k0=k[i-1]
  x[i] = x0*np.cos(2*np.pi*0.31) - y0*np.sin(2*np.pi*0.31)
  y[i] = x0* np.sin(2*np.pi*0.31)+ y0*np.cos(2*np.pi*0.31)
  z[i] = z0* np.cos(2*np.pi*0.32)- k0*np.sin(2*np.pi*0.32)
  k[i] = z0* np.sin(2*np.pi*0.32)+ k0*np.cos(2*np.pi*0.32)
  y[i]+=0.3*(x[i]**2-z[i]**2)
  k[i]+=0.3*(x[i]*z[i])
jet= plt.get_cmap('jet')
colors = iter(jet(np.linspace(0,1,30)))
plt.figure(0)
w=30
for i in range(1,30):
  fourier = np.fft.rfft(x[i:i+w])
  freqs =np.fft.rfftfreq(len(x[i:i+w]))
  color=next(colors)
  plt.plot(freqs,abs(fourier), c=color)
plt.figure(1)
fourier = np.fft.rfft(x)
freqs =np.fft.rfftfreq(len(x))
plt.plot(freqs,abs(fourier), c='k')
plt.show()


texts: why restricted isometry property constant $\delta_{2k}<\sqrt{2}-1$?

texts:
It's said that $\delta_{2k} < \sqrt{2} -1$ , the solution of the $l_{1}$ problem is that of $l_{0}$ problem.

I checked the proof of $||x^{*}-x||_{l_{2}}\leq C_{0}s^{-1/2}||x-x_{s}||_{l_{1}}+C_{1}\epsilon$

in this file but I haven't found out where the $\sqrt{2}-1$ works.

So it will be very helpful if anyone could tell me if I missed something in that paper or direct me to some relevant literature or give me some tips.


texts: Normalization of signal against reference

texts: I am trying to figure out if there is any standard way of normalizing a sequence such as a time series against a reference sequence. I am not an expert on signal processing but I was hoping someone would be able to point out a good way of doing this. Specifically, I am working on image analysis where I measure pixel intensities along linescans. This gives me two intensity profiles along space (from 2 channels). I am trying to normalize 1 channel against the 2nd one (reference).


texts: Estimate amplitude of signal between 45 Hz-60 Hz

texts:
I have to estimate the amplitude of vibration of a real world robot.

The signal is originating from a laser distance measurement, and the vibrations are somewhere between 45 and 60 Hz.

I'm using matlab.

I have 2 options to estimate the amplitude of my signal:

Apply an FFT with a Hanning window, then select the maximum
peak in the window [45;60]Hz, then calculate amplitude with this
relationship: Amplitude = FFTPeakValue/amountOfSample*2 ( I have to
multiply by 2 because signal is symmetric around
SamplingFrequency/2, and I'm looking at one side only).
Apply a pass band filter between [45;50], then calculate the amplitude with this relationship: Amplitude=maxPeaktoPeak/2;

Are both way viable according to you ? Is one of them better?


texts: What does the physical meaning of $\mathbf{U} \Sigma \ \mathbf{V}^{*} $ in the MIMO system?

texts:
SVD is a method to cancel the interference in MIMO system,we often do the SVD of channel as below

\begin{align}
  \mathbf{H} &=
  \mathbf{U} \ \Sigma \ \mathbf{V}^{T} \\
%
\end{align}

So  $ \mathbf{U}$ means the beamforming direction of receiver? and $\mathbf{V}$ means the  beamforming direction of transmitter? and  the singular value in the diagonal of $\Sigma$ means the gain of channel ?Why can we use this to decrease or cancel the noise ?


texts: Can an untilted perspective of an object be obtained from a 2D image of a tilted object using angular transformations?

texts:
I'm trying to make a height estimation algorithm that works in tandem with a bounding box object detector. The object detector can recognise the objects in images taken at pitch/yaw angular offsets.

Given these pitch and yaw offset angles are known, is it possible to apply a transformation on the object bounding box to obtain the box's equivalent dimensions in an un-tilted perspective, which would subsequently be used by a height estimation algorithm?


texts: Method for finding the best non-increasing approximation function of a given signal

texts:
I have a discrete-time signal (with noise) and I want to determine the best non-increasing function that represents it.

Of course it would be useful if the obtained function were defined in terms of a continuous variable of time, but to me this sounds too general. I think it is not a good idea to specify a particular equation with coefficients to be determined, simply because in so doing I am imposing a too restrictive space of test functions, thus possibly causing "underfitting". For example, I could impose an exponential function like $f(t) = c \, \exp(-\lambda t)$, where $c$ and $\lambda$ are parameters to be determined, but what if the given signal has another form?

For a discrete signal I think it is better try to generate a discrete approximation as well. In this case, the "fit" would be a non-increasing discrete function that is the best representation of the given general signal.
Does anybody know a method to do that?

One could use either the raw signal or then a low-pass filter firstly. As a remark I would like to say that polynomial interpolation simply does not work because of the noise.


texts: Frequency shift property of Fourier transform

texts:
Which one of the following is actually a definition of frequency shift property
$$e^{j\omega_0 t}x(t)\leftrightarrow X(j\omega - \omega_0) \tag{1}$$
$$e^{j\omega_0 t}x(t)\leftrightarrow X(j(\omega - \omega_0)) \tag{2}$$
Is frequency shift applicaple in all the situations?

Case I
$$1 \leftrightarrow 2\pi \delta(\omega)$$
if we apply frequency shift property we may obtain  $$e^{j\omega_0 t} \leftrightarrow 2\pi \delta(\omega -\omega_0)$$
which works according to result 2
Case II
$$u(t) \leftrightarrow \frac{1}{j\omega} +\pi \delta(\omega)$$
$$e^{-at} u(t)\leftrightarrow \frac{1}{a+j\omega}$$
which exactly isn't 1 or 2


texts: Explain FFT indexing with respect to bit reversal and digital reversal of fixed 5 indices

texts: Explain FFT indexing with respect to bit reversal and digital reversal of fixed indices


texts: Software to register a 2D image to a 3D volume

texts:
I have been searching for software that can handle the task of registering a 2D image to a 3D volume.  The 3D volume is a microCT scan, while the 2D image is taken from a plane passing through the sample.  I am looking for something that can locate the plane of the 2D image within this 3D volume, i.e. get the 3D translation, rotation and scale transform betwee the image and volume.

First prize is for pre-compiled free Windows software, but if that isn't possible then any libraries that can handle this task.  I see that ITK from Kitware provides registration tools, but I don't see any discussion of my specific task, which doesn't seem to a common one.


texts: Studying spectral leakage due to DFT length and frequency resolution

texts:
I guess this is a pretty basic question, but I am kinda stuck... So if I am remembering the DSP theory correctly, spectral leakage can occur when we take "inappropriate" combinations of the following:

1) FFT length

2) sampling frequency

3) fundamental frequency of input signal

Here is a MATLAB script for demonstrating this, in a case of a sinusoid input:

F0=10; %Frequency of the sinusoid
Fs=100; %Sampling Frequency
observationTime = 1; %observation time in seconds
t=0:1/Fs:observationTime-1/Fs; %time base

x=sin(2*pi*F0*t);%sampled sine wave

N1=100; %DFT length same as signal length
X1 = 1/N1*fftshift(fft(x,N1));%N-point complex DFT of x
f1=(-N1/2:1:N1/2-1)*Fs/N1; %frequencies on x-axis
stem(f1,abs(X1));

MY thought was that if the fundamental frequency is a multiple of the frequency resolution, then there will be no spectral leakage, since there is a bin that corresponds to the exact value of the fundamental frequency F0. However, it seems like I am wrong here...

Can someone please explain to me what's the flaw in my rationale? Why is there leakage in the 3rd plot? Thanks in advance


texts: Implementation of Wikipedia Equation for the DFT

texts:
I was writing a simple fourier transform implementation and looked at the DFT equation on wikipedia for reference, when I noticed that I was doing something differently, and after thinking about it felt that the wikipedia version must be wrong because it's very simple to think of a signal that when fourier transformed (with that equation) will return an incorrect spectrum: Because the equation wraps the signal around the complex plane only once (due to the $n/N$ with $0<n<N-1$), any signal that is periodic an even number of times (while wrapping the complex plane) will have no spectrum as the usual peaks (while going around the unit circle) that would appear during a DFT will cancel each other out (when an even number of them appear).

To check this I wrote some code which produced the following image, which seems to confirm what my thoughts.

"Time using equation" uses the equation
$$ X_f = \sum_{n=0}^{N-1} x_n (\cos(2\pi f t_n) - i \sin(2\pi f t_n)) $$
with $t$ a vector of time (so the time $t_n$ at which $x_n$ was sampled for example). It can be found in the function ft below.

The wikipedia equation, linked above, is copied here for reference:
$$ X_f = \sum_{n=0}^{N-1} x_n \left(\cos\left(2\pi f \frac{n}{N}\right) - i\sin\left(2\pi f \frac{n}{N}\right)\right) $$
It can be found in the function ft2.

import numpy as np
import matplotlib.pyplot as plt
plt.style.use('ggplot')

def ft(t, s, fs):
    freq_step = fs / len(s)
    freqs = np.arange(0, fs/2 + freq_step, freq_step)
    S = []
    for freq in freqs:
        real = np.sum(s * np.cos(2*np.pi*freq * t))
        compl = np.sum(- s * np.sin(2*np.pi*freq * t))
        tmpsum = (real**2 + compl**2) ** 0.5
        S.append(tmpsum)
    return S, freqs

def ft2(s, fs):  # Using wikipedia equation
    nump=len(s)
    freq_step = fs / nump
    freqs = np.arange(0, fs/2 + freq_step, freq_step)
    S = []
    for i, freq in enumerate(freqs):
        real = np.sum(s * np.cos(2*np.pi*freq * i/nump))
        compl = np.sum(- s * np.sin(2*np.pi*freq * i/nump))
        tmpsum = (real**2 + compl**2) ** 0.5
        S.append(tmpsum)
    return S, freqs

def main():
    f = 5
    fs = 100
    t = np.linspace(0, 2, 200)
    y = np.sin(2*np.pi*f*t) + np.cos(2*np.pi*f*2*t)

    fig = plt.figure()
    ax = fig.add_subplot(311)
    ax.set_title('Signal in time domain')
    ax.set_xlabel('t')
    ax.plot(t, y)

    S, freqs = ft(t, y, fs)

    ax = fig.add_subplot(312)
    ax.set_xticks(np.arange(0, freqs[-1], 2))
    ax.set_title('Time using equation')
    ax.set_xlabel('frequency')
    ax.plot(freqs, S)

    S, freqs = ft2(y, fs)
    ax = fig.add_subplot(313)
    ax.set_title('Using Wiki equation')
    ax.set_xlabel('frequency')
    ax.set_xticks(np.arange(0, freqs[-1], 2))
    ax.plot(freqs, S)

    plt.tight_layout()
    plt.show()

main()

Obviously it seems rather unlikely that I would have randomly found an error on such a high profile wiki page. But I can't see a mistake in what I've done?


texts:
Apologies if I'm using incorrect terms, and if this question's already been answered here. I have a data set which consists of a high DC signal modulated by AC signal (taken from custom-built biomedical instrumentation). Now the problem is that my signal has some low-frequency features that I'm having a hard time getting rid of using scipy. Please note that I'm a researcher but have limited practical experience in signal processing. I first do a low-pass on the signal to remove high-frequency components, and then I am trying to add a second filter to remove this low-frequency feature.

I tried doing a (butterworth order 1, cutoff 0.1Hz) high-pass filter but this removes the DC offset, which is important for me to keep (Figure 1). It also behaves in a rather odd manner.

I then tried doing a (butterworth order 1, band 0.1Hz to 0.5Hz) bandstop filter (removing low frequencies just above zero and just below  my range of interest) but this results in a strange peak in the filtered output which I am unable to explain (Figure 2 - bonus points for someone who can shed light on this).

What I want to do is level the data, and remove these artifacts, so that each peak will  be in line with the next. I've searched around but there doesn't seem to be an obvious and readily available solution to this.

To give you an idea of what I would really like to do, here is a similar treatment of the signal done in OriginPro: (Look at the blue signal - it has not had its DC offset shifted at all. How is OriginPro implementing this highpass filter? How can I replicate this behaviour in scipy?)

Why do I want to use scipy if my job's being done in Origin? Simple: (1) I want to understand the processing in greater detail so I can optimize it to not lose important parts of the signal (2) To automate the entire processing from acquisition to analysis.

Would greatly appreciate any help at all. Thanks in advance.


texts: FIR filter reverse engineering from step response

texts:
I am working on a reverse engineering problem - I am trying to reconstruct a filter.

I have:

The step response of the system. (I could calculate the impulse response by derivating)

I know:

The system consists of a low-pass FIR and a decimation unit (decimation by 4, probably).
In the step response, there is always a point in the middle of the step amplitude (the attached image)
The group delay of the entire system is around 68.

The impulse response should be symmetrical - for all obtained impulse responses (I have acquired several step responses) there are two points with similar amplitude - apart from those two points the impulse response is symmetrical.
According to the manufacturer, the filter should have 5% skirt and 104 dB attenuation.

My main task is to repeat the behaviour of the system - convolving new signals with the obtained impulse response works - still - I would like to have a model of my filter.

I would like to synthesise a filter (FIR), which similar characteristics - with group delay around 68.  (the filter should have more than decimation*68*2 taps)

I cannot obtain a filter which would guarantee this middle point. Could you recommend me a tool/approach? (currently, I'm using the MATLAB filter design tool)

EDIT. I would like to generate several filters and select the best filter. (more Bayesian approach)

Thank you in advance


texts: Why we need fourier transform of periodic signal although we have fourier series for periodic signal?

texts:
I was going through some of the basics of fourier series and fourier transform. And I came across  one topic "Fourier transform of periodic signal". But I am not able to understnad why we need to go for fourier transform because we already have fourier series to take care of periodic signal.

Any help is appreciated.


texts: How to compute transfer function (Frequency Response Function) from experimental data

texts:
I have a set of data, in the freqeuncy domain, corresponding to a structural test.

I have amplitude, phase and frequency for both the measured output and the input signal (accelerations are measured in the test). The input covers the frequency range $5-100 Hz$.

I can compute the frequency response:

RE = abs(amplitude)./sqrt(1+(tan(phase)).^2);
IM = RE.*tan(phase);
response = complex(RE,IM);

Having the frequency responses of both input and output I could compute the transfer function $H(\omega)$:

H = abs(response)./abs(g_inp); % g_inp is the frequency response of the input

The problem arises from the fact that the 2 vectors (response and g_inp) have not the same length: the first one has 849 samples and the second one 828.

How can I 'synchronize' them?

EDIT

I tried following the zero-padding method as suggested but I do not get the expected plot. Here my code:

% The input signal is a sine-sweep for which I have:
% ginz: amplitude of the input signal
% finz: frequency of the input signal

% Sine sweep in time domain
R = 4;
tz = 60/R*log2(finz/finz(1));
u_swt = sin(2*pi*((60*finz(1)/(R*log(2.))*(2.^(R/60*tz)-1))));
base_acc_t = ginz.*u_swt;

% zero-padding
acc_input = [base_acc_t; zeros(N-n,1)]; % N: length output signal (frequency domain)
                                        % n: length input signal (frequency domain)
input_freq = fft(acc_input);

H = abs(response)./abs(input_freq);

I obtain a plot which has the peak value at the right frequency, but instead of a nice curve it looks really noisy


texts: Mathematical model of equivalent-time sampling, the resulting unevenly spaced periodic signal, and its interpretation

texts:
I sample a continuous signal $s(T)$ over time. This leads to $s[t]$, which depends on several factors of which some are (pseudo-) periodic. I am interested in the effect of one such periodic factors and want to look at signal changes during this period under differing experimental conditions. I know the start and end of every period $i$ and said period slowly and slightly varies over time.

Raw signal $s[t]$

My first approach was to estimate the size of the periodic changes by looking at the power spectrum, where under certain conditions the largest peak is indeed around the frequency for the mean period which I know a priori.

Power spectrum $S[\omega]$

For each measurement I know the start and the end time of the current period $i$, say $t_{i,start}$ and $t_{i,end}$. If I calculate the "phase" or relative position in the current period $\varphi(t)=\frac{t-t_{i,start}}{t_{i,start}-t_{i,end}}$ I can reorder my data as such and plot it. Because I have many samples and the frequency of the effect is not correlated to my sampling frequency the phase increase in my reordered data averages out to pretty much linear (so the two images below are almost identical).

Reordered raw data for reordered sample number $s[n]$

Data plotted to phase $s[\varphi]$

Looking around with google I found that this is called Equivalent-Time Sampling (ETS) in some fields. However it is difficult to find a mathematical formulation for this operation, so it is difficult to see really what information I keep and lose from my data.

My question now is broad: is there a formal mathematical way to model ETS (maybe as a filter on a transformed domain?) and how can I make valid observations about the size and qualities (for example waveform) for this effect when using ETS?

Above (when calculating $\varphi$) I implicitly asume that the waveform is stretched if a period is longer, but in reality there might just be a longer delay before a certain impulse is processed again. Might it be better to use an absolute time delay since start of current period? (as $\varphi(t)=t-t_{i,start}$)

Update

So I think the data, once in equivalent time, can also be considered as unevenly spaced timeseries, since $\varphi(t)$ is unevenly spaced due to the differing period lengths. An answer here suggests methods to calculate correlation between unevenly spaced timeseries. With those methods it might be possible to make a comparison between measurements without any filtering and without assuming a certain sampling rate (seems artificially high in the reordered data). Another answer (and this related answer) suggests to use a Kalman filter to smooth the unevenly spaced time series.


texts:
Assume signal $x$, sampled at $f_s = 44100\; \mathrm{Hz}$. I tried to filter it using the Butterworth bandpass filter ( $30\; \mathrm{Hz} - 70\; \mathrm{Hz}$) of order $8$. However, as a result I get a vector with most elements being NaN (and some of them extremely small, approx.  $-2.5 \cdot 10^{306}$`).

If I try the same filter of order $6$, I get results as expected. What could be possible reason for order $8$ filter to 'explode'?

Here is the MATLAB code, just in case I made an error which I don't see:

[b, a] = butter(4, [60 / fs, 140/fs]);
x_filtered = filter(b, a, x);


texts: How to preprocess signals prior cross-correlation?

texts:
Problem Description

I have two devices which are measuring the same inputs and are using them to calculate the same (non-periodic) output signal. The process generating the inputs is deterministic.
Both devices are sampling with different rates (each > 5 kHz).
The capture duration of one device may be several seconds longer than the other one. Total capture duration is > 60s.
Since the calculation of the output signals is done with different algorithms they are looking quite different at first sight but they do look very similar after applying some sort of linear low-pass filter.

My goal is to find out the time delay between those two output signals x and y.
Proposal for Solution

Crop the longer signal such that x and y have the same capture duration.
Resample x and y such that they have the same number of samples.
Filter x and y with the same linear low-pass filter.
Remove linear trends by applying MATLAB's detrend function to x and y.
Apply MATLAB's finddelay on x and y.

Questions

May cropping the longer signal corrupt the result of the cross-correlation?

What is the best approach for resampling x and y?

As I mentioned, x and y are only looking very similar
after filtering. Is filtering signals prior cross-correlation
advisable? If yes, are there recommended filter types for that
purpose?

Is it a must to detrend data prior cross-correlation or might it be enough to remove the mean? In my case there may be multiple linear trends and it could be hard to find all breakpoints as depicted here.

Desired Time Delay between two (filtered) Output Signals


texts:
What would be the best way to take the Fourier transform of
$$
f(t)\cdot \cos\big(\pi(t-1)\big)
$$
I'm aware that when you take the Fourier Transform of $\cos(kt)$ you get two impulse at the location of $k$. And that the $(t-1)$ would create a complex exponential. But you can only apply the time-shift when both functions are time shifted (at least to my knowledge at least). Really the $\pi$ and the non-time shifted $f(t)$ are really throwing me off. Thanks!


texts: How to decide about length of vectors for estimating correlation matrices

texts:
I am optimizing a filter $f$ using mean squared error criterion.

$$MSE = f^T*(R_{xx}-R_{yx}^T*R_{yy}^{-1}*R_{yx})*f \tag{1}$$

where

$x\in \mathbb{R}^n$ input column vector
$y\in \mathbb{R}^m$ output column vector and $y=Hx+AWGN$. $H\in \mathbb{R}^{m\times n}$ is ISI channel response convolution matrix.
$R_{xx}=E[x^{}x^T]$,  $R_{yx}=E[y^{}x^T]$, $R_{yy}=E[y^{}y^T]$

I know that $f$ can be found by imposing some meaningful constraints on (1).

My question is about finding the right $n$ for estimating $R_{xx}, R_{yx}, R_{yy}$. I know that as I increase $n$ (the dimension of my input vector $x$), these estimates improve. Therefore, if I fix $m$ (dimension of $y$) and plot $MSE$ on (y-axis) vs $n$ on (x-axis), I have a decreasing curve from left to right.

What I don't know is how to find the right value of $n$ when memory of ISI channel is $\mu$.


texts:
Hyper-spectral images are images with many bands $(\ge 100)$ while the normal color images have only three bands $\rm (r,g,b)$.  It is not difficult to simulate an artificial gray-scale image with the following Matlab codes:

width = 200;
height = 200;
img = ones(height,width)*10;
img(floor(height/4):floor(height*3/4),floor(width/4):floor(width*3/4))= 200;
noise = 5+10*randn(height,width);
img = img+noise;
img(img<0)=0;
img(img>255)=255;
figure; imshow(img,[]);

The following simulation image can be obtained:

In the simulation image, the object is a square in the center and the rest is the background. The object is of gray-value 200, and both object and background are contaminated by Gaussian noise. It is also possible to perform blurring before adding noise. The blurring is used to simulate the system point spread function.

Then comes to my question: how can we simulate a hyper-spectral image?

For the time being, my solution is as follows:

Find a hyper-spectral profile from the hyperspectral database for a certain object (this kind of database exists?)
For each band get the object's gray-value from the profile, and treat this band as if it is gray-scale image when creating simulation image.

I am not sure whether this is the right way of generating simulation images. Any ideas?


texts: Minimum time delay that can be estimated between two sensors

texts:
Consider a 1-dimensional toy problem. I have two sensors at different points along the $x$-axis. Somewhere away from the sensors, a disturbance is created which travels towards the sensors, first reaching sensor 1 and then sensor 2 at a later time. Both sensors produce samples over the entire time window, $W$. Both sensors are sampling at the same rate $f$. Using cross-correlation or other methods, what is the smallest time delay I can resolve?

My hunch is that the answer is given by the Cramer-Rao bound. If so, can someone give a simple intuition on how this bound is derived?

Thanks!


texts: Relation between FIR filter's transition width and phase delay

texts:
How does the transition width of a FIR filter relate to phase delay at the output of the filter?

I have been trying to find information on the subject for days.


texts: Echo Removal from PCM 16 bit arrays (time frames)?

texts:
I am keen to know the basics to remove echo knowing FarEnd Signal and NearEnd Signal. I have read research papers in which its commonly said:
y(n) = (hT)x
where h is the room impulse response and T is transpose and echo signal estimated is y(n).
Now, what I really want to understand is the computation of room impulse response. Given two-time frame arrays of PCM 16 bit. Let's say:
Near_end_sig = [x1,x2,x3,x4,....,xN] (one time frame)
Far_end_sig = [y1,y2,y3,....,yN] (one time frame)
Here I want to know in this case working on digital acoustic signals (PCM-16 bit arrays of size let's assume 1024). How do I compute room impulse response h and how do I estimate echo signal? My main concern is the computation of room impulse response and the echo signal estimation based on the time frames.
Regards,
Khubaib Ahmad


texts: What does it mean that the DFT equals the Complex Fourier Coefficients for even frequencies?

texts:
I have a periodic signal, with period $1$

$$x(t) = \begin{cases}
1 \qquad & 0 \le t - \lfloor t \rfloor < \tfrac12 \\
0 \qquad & \tfrac12 \le t - \lfloor t \rfloor < 1 \\
\end{cases}$$

$\lfloor t \rfloor = \operatorname{floor}(t)$ is the floor() function, returning the largest integer no greater than the argument $t$.

$$ x(t+1) = x(t) \qquad \forall t \in \mathbb{R} $$

The complex Fourier series for $x(t)$ is

$$ x(t) = \sum\limits_{k=-\infty}^{\infty} c_k \ e^{i 2 \pi k t} $$

The complex Fourier coefficients are

$$\begin{align}
c_k &= \int_{-1/2}^{1/2} x(t) \ e^{-i2\pi k t} \ \mathrm{d}t \qquad \qquad k \in \mathbb{Z} \\
    &= \int_{0}^{1/2} 1 \ e^{-i2\pi k t} \ \mathrm{d}t \\
    &= \tfrac{1}{-i2\pi k } \big( e^{-i\pi k} - 1 \big) \\
    &= \tfrac{i}{2\pi k } \big( (-1)^k - 1 \big) \\
\end{align}$$

and are $0$ for even $k$.

Now, we sample $x(t)$ at $N\in 2\mathbb{N}$ time values,

$$\begin{align}
x[n] &= x(t_n) \\
     &= x\left(\tfrac{1}{N}n\right)
\end{align}$$

where $\tfrac{1}{N}$ is the sampling period, $N$ is the sampling frequency, and $t_n = \frac{n}{N}$, with $n=0, N-1$.

Let $X[k]$ denote the DFT of this finite sequence $x[n]$.

$$ X[k] = \sum\limits_{n=0}^{N-1} x[n] \ e^{-i2\pi nk/N} $$

One can show that $X[k] = 0$ for even $k$ as well, so $X[2k] = c_{2k}$, but $X[2k+1] \neq c_{2k+1}$.

I'm trying to come up with an intuitive explanation as to why $X[2k] = c_{2k}$, but $X[2k+1]\neq c_{2k+1}$. Obviously, we shouldn't really expect them to be equal in general since $X[k]$ is really just a Riemann sum approximation to $c_k$ with $N$ intervals, but in this case it seems there may be an explanation, since $X[2k] = c_{2k}$.

So far, all I can come up with is that since $x(t)$ isn't bandlimited, the DFT of the discrete sample of $x(t)$ is essentially trying to "fit" a bandlimited signal to the samples $x[n]$, and so for some reason this doesn't have any content at frequencies $2\pi (2k) = 4\pi k$ for any $k\in\mathbb{N}$.

Is there any specific reason for this? From the above argument, I have a feeling it has to do with aliasing, but I can't exactly make the connection.


texts:
Fourier series in continuous time domain while representing $a_k$ in rectangular form   $$ a_k = B_k + jC_k$$
But when using the value of $a_k$ in the main equation: $$ x(t) = a_0 + 2\sum^{+\infty}_{k\ =\ 1} [B_k\cos k\omega_0t - jC_k\sin k\omega_0t]$$

I want to ask from where $-$(minus) sign comes as there was no minus sign
main Equation
$$ x(t) = a_0 + 2\sum^{+\infty}_{k\ =\ 1} 2\mathrm{Re}\left[ e^j(k \omega_0t +\theta_k)\right]$$


texts: Filter design with 8-tap latency on recursion

texts:
I'm curious about the feasibility of designing a noise shaping filter (but might generalise to any recursive filter) with the constraint that the most recent output samples aren't available for several iterations of the filter.

The use case I have in mind is reducing the word length of an audio stream sampled at 48kHz, and the filter would try to keep most of the quantisation noise out of the 200Hz-8kHz band (where human hearing is more sensitive); but with the additional constraint that the most recent seven outputs are not available to the filter.

I'm applying this constraint so that the output can be produced several samples at a time using a mixture of SIMD and instruction-level parallelism -- working with a notional four-lane SIMD and two-fold unrolling.

There's a lot of work for each sample which could be performed in parallel (scaling, quantisation, dither) if not for the recursive factor.  Normally one might hope to use SIMD for the multiply-accumulate of the coefficients, but that's not a great proportion of the work if the filter has only a handful of taps.

So I'm wondering if it's actually possible to design a filter that doesn't get in the way, here, but is still effective for its intended purpose.

Stereo would obviously offer an opportunity to halve the delay, but avoiding that assumption is preferable, as it makes mono a failure case and multi-channel becomes complicated.


texts: What are the units of my data after an FFT?

texts: Magnetometer measures the derivative of the magnetic field, or dB/dt, with an  output in microvolts (mV). The Sampling rate is 128 Hz, so if we collect data for 2 minutes, $2 \times 60 \times  128=15360$ points (discrete case). When I perform an FFT on these time series, what will the units of amplitude density be  after the transform?


texts:
Let's say I have an arbitrary signal which has the following frequency spectrum and a bandwidth of BNF:

Now I use Amplitude Modulation to transport the signal. The frequency spectrum will change accordingly:

Now I know that the required bandwidth doubles. I don't understand where the doubled bandwidth applies as I can't transport negative frequencies.


texts: Find state space model from transfer function

texts:
Let's suppose we have:

G(s) = (s+1)/(s^2-2s+1)

how can we find the state space representation of the transfer function:

x_dot = x2

x2_dot = 2*x2-x1+u

where u is an arbitrary input.

I am very new to this topic, so a detailed answere would be great ! :)


texts: Phase tracking in communication receivers

texts: Why does Costa's loop PLL bandwidth need to be narrow? Text books mention that and yet give no explanation


texts: An Interesting Model with Unknown Orthogonal Design Matrix

texts:
Suppose the   multivariate one-way anova  model for the raw data , i.e.
$$
\label{Example_model_1}
\mathbf{y}_{ij}=\mathbf{\mu}+\mathbf{z}_i+\mathbf{e}_{ij}, ~~ i=1,\ldots,m,~~j=1,\ldots,n_i,~~~~~~~~~~~~~~~~(1)
$$
where  $\mathbf{\mu}$ is the unknown overall mean of $\mathbf{y}_{ij}$,
$\mathbf{z}_i$ is the unknown random effect at the $ith$ level,  $\mathbf{e}_{ij}$ are the errors,  and
$$
\mathbf{e}_{ij}\mid \mathbf{\Lambda}_0 \overset{iid}{\sim} N_k(\mathbf{0},\mathbf{\Lambda}_0), ~~
\mathbf{z}_i \mid \mathbf{\Lambda}_1 \overset{iid}{\sim}  N_k(\mathbf{0},\mathbf{\Lambda}_1).
$$
Howerever, the raw data is not observed; instead $\mathbf{y}_{ij}^* = \mathbf{\Gamma} \mathbf{y}$ is observed, where $\mathbf{\Gamma}$
is an unknown orthogonal matrix (i.e., the vector observations have been randomly rotated, with
a common unknown rotation).
So the model (1) becomes
$$
\label{model_1_1}
\mathbf{y}^{*}_{ij}=\mathbf{\mu}^*+\mathbf{z}_i^{*}+\mathbf{e}_{ij}^{*}, ~~ i=1,\ldots,m,~~j=1,\ldots,n_i,~~~~~~~~~~~~~~(2)
$$
where $\mathbf{\mu}^*=\mathbf{\Gamma}\mathbf{\mu},$ $\mathbf{e}_{ij}^*\mid (\mathbf{\Lambda}_0,\mathbf{\Gamma})
\sim N_k(\mathbf{0},\mathbf{\Gamma}\mathbf{\Lambda}_0\mathbf{\Gamma}'),$ and
$\mathbf{z}_i^* \mid (\mathbf{\Lambda}_1,\mathbf{\Gamma}) \sim N_k(\mathbf{0},\mathbf{\Gamma}\mathbf{\Lambda}_1\mathbf{\Gamma}'). $

Our aim is to estimate $(\mathbf{\mu},\mathbf{\Gamma},\mathbf{\Lambda}_0,\mathbf{\Lambda}_1).$

Question:
I can't find the relative background of the above model. Please give me some reference or suggest.


texts: Changes of phase in high frequency makes the speaker "clicks"

texts:
I modulate some data over sound in a 17.7Khz and 17.8Khz. It works like an FSK with a pulse duration of 0.015sec. When playing a string the speaker gives a very strange noise that it seems to be the phase 180 degrees change.

Is there any way to reduce this noice ?


texts: Normalized Frequency to rads/sec, Hz and discrete/digital frequency

texts:
I have a signal sampled at 16000 Hz
I need to convert the normalized frequency of Matlab 0.325 to radians/sec Hz and Discrete frequency.
Can anyone explain this to me?
And btw Matlab plots the normalized frequency as $\pi$ rads/sample. So we have to divide the normalized frequency with $\pi$?


texts:
Both Correlation and Convolution are displacement function, i.e they are used to slide the filter mask across the image.

Convolution is same as correlation except that the filter mask is rotated 180 degree before computing the sum of products.

We mostly use convolution. why it is widely used? why convolution is preferred over correlation?

We have,

 Convolution in time domain   = Multiplication in Frequency domain

Do we have any similar relation for correlation ?

Correlation in time domain  =  _______________ in Frequency domain ?


texts: Octave demodmap function - simple assist

texts:
I'm a beginner with Octave.  I'm trying to de-modulate a simple FSK file (wav file of alternating 2KHz and 4KHz tones at 500 bits per second, into an array of 1s and 0s. I read the documentation for function demodmap in the communications package, but I must be understanding something wrong. I need a simple example that works.

data = audioread("monoshort8008n9.wav");
z = demodmap( data, 500, 44100, "fsk", 2);

Two errors occur. First 44100 is not a multiple of 500 as required by the function, so I tried 44000 instead (knowing it is incorrect). Then I just get a z array of zeroes.

Questions: Most of the time sample rate is not a perfect multiple of the symbol rate. Is this really a requirement of Octave?  If so, how do I "stretch" or "resample" my data at 44000 instead of 44100, which is what Audacity really did?

But even then, the function is returning all zeroes.  Any ideas?


texts: FFT of white noise with different sampling rate

texts:
I have a question on how I should interpret the white noise power level (noise floor) obtain from FFT for different hardware sampling rates. I realized if I sample the same noise at different rates (100 MHz, 10 MHz), even though I have down-converted (FIRfilter and decimate) them to the same, lower rate (1 MHz), I will obtain a different noise level if I FFT the data.

For example, with 1 sec worth of data, I obtain 100 million points sampled at 100 MHz. By down-converting it to 1 MHz and doing FFT, I obtain a noise floor of -80 dBm. If now I reduce my hardware sampling rate to 10 MHz and down-convert it to 1 MHz, the FFT result will give me -70 dBm of noise floor which is higher than when I sampled at 100 MHz.

May I know what the cause is?

My initial guess was the broad-band noise caused aliasing and they added up. Does it mean I should always filter out physically the noise above double my sampling frequency?


texts: Computations in Grayscale morphology dilation

texts:
I red this link, but unfortunately i still don't understand the process...

By the way, Dilation process is performed by laying the structuring element H on the image I and sliding it across the image in a manner similar to convolution.

The grayscale morphological dilation formula is written as follows :

$[I ⊕ H](u,v) = \displaystyle\max_{(i,j)∈H}\{I(u-i,v-j)+H(i,j)\}$

If we assume a greyscale image I of size 2x3, such that I=[1 2 3; 3 7 2]; and a structuring element H of size 1x4, such that H= [0 0 7 7].

What are u,v and i,j? Please can someone explain to me in details how the formula is computed at each element of I in order to obtain the new image dilated_I ? Because i still don't understand well the concept..

What was meant by a structuring element that contains only the origin ? i didn't understand well the concept behind grayscale dilation.
Any help will be very appreciated.


texts: Sphere decoding for non-binary modulation

texts:
I have a simple question on how to implement SD for non-binary modulation, let say QPSK. Assume we transmit a complex symbol $\overline{a}$ drawn from a QPSK constellation. On a vector form, we can write the receive signal as
$$\overline{y} = \overline{G} \: \overline{a} + \overline{w},$$
where $\overline{G}$ is the channel matrix and $\overline{w}$ is the complex noise samples. To detect $\overline{a}$ from $\overline{y}$, we can use maximum likelihood (ML) detection or to simplify the detection we use sphere decoding (SD).

To use SD, the above equation can be rewritten in a real-valued form as
$${y} = {G} \: {a} + {w},$$
where $y = \begin{bmatrix}
\Re\{{\overline{y}}\}\\
\Im\{{\overline{y}}\}
\end{bmatrix},$ $G = \begin{bmatrix}
\Re\{{\overline{G}}\} & - \Im\{{\overline{G}}\}\\
\Im\{{\overline{G}}\} & \Re\{{\overline{G}}\}
\end{bmatrix},$ $a = \begin{bmatrix}
\Re\{{\overline{a}}\}\\
\Im\{{\overline{a}}\}
\end{bmatrix},$ $w = \begin{bmatrix}
\Re\{{\overline{w}}\}\\
\Im\{{\overline{w}}\}
\end{bmatrix}$, $\Re$ and $\Im$ are the real and imaginary parts of complex numbers. Now, the SD can be applied directly with no problems no the new problem of double the dimension.

My question is can we apply two parallel versions of the SD algorithm, one to detect the real part and the other to detect the imaginary part?


texts: continuous chain of ones or zeros for transmit

texts:
If my data has long chain of 1's or 0's. Can the data be modulated and transmitted as it is. Will this have any effect of  the receiver adc as this may increase dc value of the signal or any other block in the receiver chain?
should the data be scrambled?


texts: Negative instantaneous frequency with hilbert transform using scipy hilbert

texts:
When I produce a chirp signal with a duration of 2 sec, sampling frequency of 200Hz, with frequencies ranging from 20Hz-100Hz, I get negative instantaneous frequencies using the scipy(Python) Hilbert transform. Here is my signal code:
    self.duration = 2
    self.fs = 200.0
    samples = int(self.fs * self.duration)
    self.time = np.arange(samples) / self.fs
    carrier = chirp(self.time, 20, self.time[-1], 100)
    self.signal = carrier * (
        1 + 0.5 * np.cos(10 * 2.0 * np.pi * self.time)
    )

Here is my Hilbert transform code:
demean = data - np.mean(data)
hil = hilbert(demean)
inst_amplitude = np.abs(hil)
inst_phase = np.unwrap(np.angle(hil))
inst_frequency = np.diff(inst_phase) / (2 * np.pi) * sampling_frequency
inst_frequency = np.insert(inst_frequency, 0, 0)
regen_carrier = np.cos(inst_phase)

Here is the plot of the Hilbert properties:

Why are there negative instantaneous frequencies? I have been Googling the issue with zero luck finding an answer. Any input would be greatly appreciated.
UPDATE: I applied the exact same Hilbert transform code as above to a "real life" dataset with a sampling frequency of 800Hz. I am getting negative instantaneous frequencies in the middle of the data, not just at the edges. I don't think I am allowed to show the data as it is proprietary, so sorry on that front. Any insight as to why I am seeing negative frequencies in the middle of my data and not just at the edges?
UPDATE #2 I went ahead and generated a synthetic chirp signal with multiple frequency sweeps using scipy.signal.sweep_poly and I am getting negative instantaneous frequencies in the middle of my data.
Here is the code for my new signal (duration of 50sec and sampling frequency of 100):
    self.duration = 50
    self.fs = 100.0
    samples = int(self.fs * self.duration)
    self.time = np.arange(samples) / self.fs
    p = np.poly1d([0.025, -0.36, 1.25, 2.0])
    self.signal = sweep_poly(self.time, p)

Here is the plot of the new signal's Hilbert properties:

Here's the same image, but zoomed in where there are negative instantaneous frequencies in the middle of my data:

Why are there negative instantaneous frequencies in the middle of my data? There are definitely no edges in the middle of my data, so it's not an edge effect.
UPDATE #3 I added random noise to the exact same signal using the following code:
    self.duration = 50
    self.fs = 100.0
    samples = int(self.fs * self.duration)
    self.time = np.arange(samples) / self.fs
    p = np.poly1d([0.025, -0.36, 1.25, 2.0])
    carrier = sweep_poly(self.time, p)  # Carrier signal
    noise = np.random.normal(0, 1, carrier.shape)
    self.signal = carrier + noise

Here is the plot of the noisy Hilbert transform properties:

Even more negative frequencies! Why is this happening?
UPDATE #4 I decided to multiply 2 sine functions with two different constant frequencies together to see what I get. There is no noise in this instance. Here is my code:
    self.duration = 2
    self.fs = 100
    samples = int(self.fs * self.duration)
    self.time = np.arange(samples) / self.fs
    self.signal1 = np.sin(1 * 2.0 * np.pi * self.time)
    self.signal2 = np.sin(9 * 2.0 * np.pi * self.time)
    self.signal = self.signal1 * self.signal2

Here is the plot of the resulting Hilbert properties:

If I change the sampling frequency, the magnitude and sign of the instantaneous frequency spikes change! Here are a few examples (the sampling frequency used is in the title of the image):

Something about the phase calculation is unstable and the sampling frequency seems to alter things. No I am even more curious as to what is happening.
Thanks!


texts:
I hope someone can help me understand the meaning of variance in this case. I am reading an old research paper by Ephraim and Malah (for my DSP project) and it says:

Knowing that $D_k$ is the FFT of white noise, what I am trying to find out is $\lambda_d(k)$. Can anyone explain to me the meaning of $E\lbrace|D_k|^2\rbrace$? I think it is the variance of a particular $D_k$ with respect to the mean, but since this is white noise, shouldn't it be zero?

Please help, I am stuck!!!!

Thanks.


texts:
From Modern Control Engineering 5th Edition page 469

The transfer function of a plant controller by a proportional controller is given by

$$G(s) = \frac{K(s+0.5)}{s^3 + s^2 + 1}$$

In the book $G(j\omega) / K$ is plotted and from that they infer that for $ 2 < K$ the system is stable.

I plotted the nyquist plot of the plant itself

and tried getting the answer consistent with the procedure from our lecture.

I know that I need 2 counter clockwise encirclements of $(-1,0j)$ so I want to see for which gain the plot is shifted past -1.

$$\frac{-1}{K} > -0.5$$

$$K > 2$$

So far so good. They did not compute the upper limit of $K$.

I think I can see that the encirclements occur between -0.5 and 0 from the plot.

So next I need to calculate

$$\frac{-1}{K} < 0$$

$$\frac{1}{K} > 0$$

$$K > \infty$$

and this can't be. I assumed to get to $$K < \infty$$

I think my main problem is that I don't exactly know when to use < and > and which points (-0.5,0,0.5) I should choose for my calculations. Also how would I proceed if there was only 1 RHP and I'd need only 1 ccw encirclement.

Could someone clarify the general procedure for this or guide me to a ressource that explains it for beginners?


texts: Maximum phase with the same magnitude response

texts:
I'm trying to solve an exercise in a general way and I can't find if my answer is correct.

Given this frequency response:

$$H(z) = 4 + 2\sqrt2 z^{-1} + z^{-2}$$

I need to find a frequency response that has maximum-phase response $H_{1}(z)$ such that $\lvert H(z)\rvert = \lvert H_{1}(z)\rvert$.

First of all, I write:

$$H(z) = \frac{4z^2+2\sqrt2 z + 1}{z^2}$$

So I see that I have zeros in:

\begin{align}
z_{01} &= \frac{\sqrt2}{4}+j\frac{\sqrt2}{4}\\
z_{02} &= \frac{\sqrt2}{4}-j\frac{\sqrt2}{4}
\end{align}

And poles in:

$$z_x = 0 \text{ (double)}$$

Now, I know that for a maximum phase system I need that all the poles and zeros have to be outside the unit circle, but from Oppenheim-Schafer I also know that I have to add a scalar factor for the minimum-phase case, so I suppose I also have to add that scalar factor now.

So I say that this:
$$H_1(z) = k \left(z- \frac{1}{z_{01}}\right) \left(z- \frac{1}{z_{02}}\right)$$
(I used that I have both poles in $z=\infty$, is that right?)
has to be the one that makes $\lvert H(z)\rvert = \lvert H_{1}(z)\rvert$.

But I can't find the value of $k$ (scalar factor), because when I write $z=e^{j\omega}$ and use the magnitude definition (magnitude = square root of the squared real part plus the squared imaginary part) I find that $k=1$ and $k=-1$.


texts:
I'm trying to code the HPS (harmonic product spectrum) algorithm and the problem I'm facing is that using my guitar, the fundemental frequency (Here it's 82Hz) just isn't there, so is the 5th Harmonic at 410Hz (Actually it's the lowest local point).
In my algorithm, I did the step of finding the max of all (Mul(harmonic-multiplications))
What pre/post-process should I do to get to 82Hz here.


texts: Time scaling of discrete-time sequences and the DTFT

texts:
In the second edition of Signals and Systems by Alan Oppenheim, he discusses the DTFT of a "time-expanded" sequence that is effectively a slowed down version of the original sequence and can be thought of as the original sequence $x[n]$ with $k-1$ zeros placed between each sample. However, the DTFT for the sequence $x[an]$ where $a$ is an integer is not discussed. Why is that? Is it because information is lost when we talk about $x[an]$ for example (since it only contains samples evaluated at multiples of $a$)?
One other thing that I noticed is that unlike the CTFT, Oppenheim only discusses differentiation in the frequency domain, unlike the CTFT where both integration and differentition in the frequency domain were presented. Why is that?


texts:
I'm learning of digital filters and found a question asking why a filter with a impulse response

$$h[n]= \begin{cases}
      \frac{1}{n^2+3n+2} & n\geq 0 \\
     0 & n<0 \\
   \end{cases}
$$

cannot be implemented. My only guess is because of the infinite amount of computation required, but the problem in the book states it is not only beacuse of this. Also, which would be good alternative filters with a finite complexity such that thay would have a close frequency response to this filter?

Cheers


texts: How to show that this filter is a HP filter

texts:
i got a system with transfer function given by:

$$H(\omega)=1-e^{-j\omega}$$
I already plot it, and that seems to be a periodic function with $H(0)=0$, $H(\pi)=2$, , is that enough to show that this is a FIR bandpass filter. thank you very much


texts:
I have a line scanner output that is scanning the geometry of a surface. Typically the output looks something like this:

Where the x and y axis refers to the coordinate system in mm. (ignore the red dots). This data is essentially a 2d point cloud.

So what i am trying to achieve is to detect the corner points of this scan. I try to start using image processing techniques, for example using a harris corner detector.

However, the issue with this is that this data has a resolution of up to 0.01um vs a possible area size of 140mm X 100mm of all possible data points. If i were to convert this to 2d image, it would be a very large file, with most of the pixels being empty.

My question is that is there a better way to go about processing this type of data; or is there an alternative way to apply techniques such as the harris corner without resorting to conversion into an image.


texts: Implementation of a Variable Fractional Delay with Lagrange Interpolation using Farrow Structure

texts:
I am writing a C++ simulation software working in time domain. I generate regularly sampled data, and need them to be delayed, in "real-time", by a variable fractional delay.

This is a pretty common problem, and very close to all the "resampling" issues. After a short bibliographic search, it seems that the most commonly used method is Lagrange Interpolation with a Farrow structure. Based on J. O. Smith, "Interpolated Delay Lines, Ideal Bandlimited Interpolation, and Fractional Delay Filter Design", I used the following implementation:

which includes cascading filters with fixed coefficients $1+z^{-1}$, modulated by cascading factors that depend on the variable delay. This can be seen in the equation

$$
\hat{H}_\Delta(q^{-1}) = \sum_{n=0}^N \hat{h}_\Delta[n](q-1)^n.
$$

My implementation in C++ works well for interpolation orders up to 51.
Unfortunately, for larger interpolation orders the numerical errors increase rapidly and become dominant. After a quick inspection, it seems that the buffered samples, in the $1+z^{-1}$ cells of increasing index, have an increasing order of magnitude. It helps when I use long double, but this becomes insufficient when the order is of the order of 60: we loose a lot of numerical precision when we combine floating-point values of around $10^{10}$ with values of the order of $1$.

Do you know a variation of the Farrow structure or another implementation for Lagrange interpolation that works well for large orders, i.e. that remains computationally efficient while requiring only reasonable numerical precision in the memory buffers?

Thank you in advance for your time :)


texts: What exactly does compression say about correlation of data?

texts:
I've been using the following formula on various empirical data $d$, to obtain a correlation factor $c_f$:-

$$ c_f = { |C(d_s)| \over |C(d)|} $$

where $C$ is a compression function like bz2 or zip, and $d_s$ is the data having been randomly shuffled. The purpose of the random shuffle it to destroy any byte wise auto correlation that might exist in $d$.

A typical value for $c_f$ might be 2.7.  What might this indicate viz a viz any possible correlation within the data?  Is it at all possible to convert 2.7 to a lag figure of $n$ bits, even loosely?

Example figures are:-

$|d|$ = 921,600 bytes

$|C(d)|$ = 201,000 bytes

$|C(d_s)|$ = 543,000 bytes


texts: Number of zeros of a sum of Shah functions by applying Rice's formula?

texts:
There is a Dirac pulse train following the scheme of the Shah function (or $\delta$-cumb function) with its Fourier series of the form:

$$\varsigma(t,T)=\sum_{n=-\infty}^{\infty}\delta (t-nT)=\frac{1}{T}\sum_{n=-\infty}^{\infty}\exp\left(\frac{i 2\pi n t}{T}\right)$$

here $n\in \Bbb Z$; $t$ time and $T(\in \Bbb N)$ the period between the pulses, which occur exactly at the levels of the integers.

The sum of $\varsigma(t,T)$ over $T$ in a given range $[1,T_{\max}]$ and $n$ in a range $[1,T]$ (instead) is accordingly:

$$S(t,T_{\max})=\sum_{T=1}^{T_{\max}}\varsigma(t,T)=\sum_{T=1}^{T_{\max}}\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t}{T}\right)$$

And the zeros (at the level of some integers) of this function with respect to $t$, hence:

$$\sum_{T=1}^{T_{\max}}\frac{1}{T}\sum_{n=1}^{T}\exp\left(\frac{i 2\pi n t_0}{T}\right)=0$$

The Question is to identify the number of such zeros $\mathcal Z$ in a range of $t\in[T_{\max},T_{\max}^2]$ by applying Rice's formula, while assuming that the frequencies $1/T$ are rationally independent?


texts: What distribution is the easiest to compress?

texts:
I'm currently playing around with some compression algorithms and I'm asking myself if there is a type of data distribution / noise distribution that is easier to target with quantization (meaning less distortions at same rate). To my understanding i.d.d. Gaussian is the upper bound on "compression difficulty".

Are there distributions that especially efficient to compress, maybe via quantizers designed for those (I know about LLoyd-Max but I'm looking for something that works generally well with a known distribution and doesn't need to be optimized samplewise)? Does it make sense to transform the data into such a distribution prior to quantizing it?


texts:
What are the DAC and ADC architectures in SDRs for inphase and quadrature phase channels?

Do we have separate ADC and DAC for each channel or there is some clever way that they are  multiplexed to reduce area.

Usually in diagrams one DAC and one ADC is shown. But then I wonder how the inphase and quadrature phase information will be separated. Would appreciate if we can get a architectural diagram with an explanation


texts:
Any function $f \in PW_{\pi}$ (where $PW_{\pi}$ is the Paley-Wiener space), can be expanded in terms of the orthonormal basis $\{e^{i n x}\}_{n\in \mathbb Z}$ as
$$\hat{f}(x)=\sum_{n\in \mathbb Z}\left\langle\hat{f},e^{i n x}\right\rangle_{L^2(-\pi,\pi)} e^{i n x}$$
where $\langle g, h\rangle_{L^2(-\pi,\pi)}=\frac{1}{2\pi}\int_{-\pi}^\pi g(x) \overline{h(x)} dx$. Taking the inverse Fourier transform in previous equation, we obtain the Whittaker-Kotelnikov-Shannon (WKS) sampling theorem,
$$f(x)=\sum_{n\in \mathbb Z}f(n) \operatorname{sinc}(x-n), \ \ \ x\in \mathbb R$$
where $\operatorname{sinc}(x)$ is the normalized sinc function. My question is the following. Sampling theorem, at least in the version that I wrote above, is related to $L^2$ space and so, to scalar product:
$$\langle g, h\rangle_{L^2(-\pi,\pi)}=\frac{1}{2\pi}\int_{-\pi}^\pi g(x) \overline{h(x)} dx$$
Instead, usually, for power signals, we define the inner product to be
$$\left\langle g\,,\,h\right\rangle = \lim_{T\rightarrow \infty} \frac 1 {2T} \int_{-T}^T g(x)\overline{h(x)}\,dx \tag{1}$$
Is there a version of Sampling theorem related to scalar product (1)? Is it used in signal processing?

References and answers are welcome.

(I premise that I'm not an expert on the signal theory, therefore I do apologise if this question isn't much precise.)


texts: I am looking for a implementation of Hilbert-Huang Transformation. Ideally, in Java, but any language with similar syntax would prove useful. Could someone point me to a specific implementation, preferably in Java?


texts: What happens when linear phase is critical?

texts: I am studying the topic digital filters. And in the slides it said,"if linear phase is not critical, Infinite Impulse Response (IIR) filters yield a much smaller filter order for a given application". So what does it mean ?


texts: Why is anisotropic diffusion useful in image processing?

texts:
I am working on anisotropic diffusion and the two coefficients proposed by Perona & Malik.

I wanted to know what is the use of diffusion in image processing? Why is anisotropic diffusion important and in which fields is it commonly used?


texts: Number of taps needed in an FIR filter to remove DC

texts:
An FIR filter is being used to remove the DC of an ECG. The sampling rate is 500Hz. How many taps would said FIR filter require (theoretically and practically) to filter out the DC.

P.S. the filter would be a high-pass filter and would not remove the 1Hz fundamental frequency component of the ECG.


texts: Why is sum of values of an edge detection filter zero?

texts:
Usually the filters for detecting edges have the sum of all the values in the filter matrix equal to 0.
For example: [-1 -2 -1, 0 0 0, 1 2 1].

Why is it so ?


texts:
Prove that $1+2\cos(\omega)$ is the same as $\frac{\sin(\frac{3\omega}{2})}{\sin(\frac{\omega}{2})}$

Knowing that
$$p[n]=\sum_{k=-1}^1 {\delta[n-k]}=\delta[n-1]+\delta[n]+\delta[n+1]$$

Doing the FT yields
$$\begin{align} P(e^{j\omega}) &=\sum_{k=-\infty}^\infty p[n]e^{-j\omega n} \\
&=e^{-j\omega}+1+e^{j\omega} \\&=1+2\cos(\omega) \end{align}$$

But, how does this fit into the general formula of $\frac{\sin(\omega(M+\frac{1}{2}))}{sin(\frac{\omega}{2})}$ where, in this case my $M=1$?

Thank you.


texts: How to compute autocorrelation with exponential IIR filter efficiently

texts:
I have a discrete dataset (called V). I want to compute the autocorrelation of V at multiple time horizons. Normally I know that I can use the IFFT of the PSD. But in my case I need it with exponential kernels.
Let X_s = V * e^(-st) for t >= 0. (Convolute V with exponentially decaying causal filter).
I want to compute <V, X_s> for many s. Each s is a different mix of every frequency so I don’t see how to do this faster than len(V) * num_s time, even with the Fourier transform.
How does one do this more efficiently the way you would with a non-exponential autocorrelation function? In the real world my time series is multivariate and with dimension 3000 and shape(V) = 1,000,000 x 3000. So my output should have num_s * 3000 * 3000 size.
Thanks!
James


texts: What is the scaling function and wavelet function at wavelet analysis?

texts:
I'm trying to looking the meaning and functionality about scaling function and wavelet function at wavelet analysis. I have googling already. But I can't find and understand the meaning.

What does those affect to analysis?
Does anyone help me to understand roughly?


texts:
I have one sample for a signal. This sample is a vector of length 384. I need to calculate the correlation matrix for this signal,So I need many samples for the same signal.
How can i generate these samples form the given sample using Matlab?


texts: Comparing two signal outputs from two clients on a networked computer game to identify lag

texts:
Basically, I'm running a server and two clients of the Unreal Engine and trying to figure out if there is lag in player movement on one of the clients (the clients are different in  way that's not important). For each client I am sampling the (X,Y,Z) location in world space every frame. Here's a sample from one player of frames 577-580:

[0037.70] ScriptLog: @@ player  577  Loc  -262.53,565.64,79.15
[0037.71] ScriptLog: @@ player  578  Loc  -261.30,563.43,79.15
[0037.72] ScriptLog: @@ player  579  Loc  -260.14,561.10,79.15
[0037.73] ScriptLog: @@ player  580  Loc  -259.15,558.78,79.15

X,Y,Z are independent. So, I have this data for two players and one player does not have lag, and so serves as the reference signal. However the second player is visibly laggy. I think a good informal definition of lag would be "large values between consecutive frames", but how do I make this notion a little bit more rigorous and also testable in a signal sense?

Here's an image of the x coordinate of the two signals. The blue player is the normal player, and the red player is the laggy player.

Edit: I'm not looking for the time delay between the two signals. Player 2 is programmed to follow Player 1 only to ensure that the trends are the same. If you zoom into the graph above, you can see that the blue line is smooth and the red line is jagged. It's the "relative jaggedness" of the red line that I want to quantify and test.

I should also mention the goal is to identify lag so that the game can start running alternative emergency routines to mitigate the lag.


texts:
I'm trying to calculate the SNR of a satellite signal using GRC, but all of the 4 offered estimators either have a lot of -NaN values or are ~35dB throughout the entire recording.
(I just connected the file source to the MPSK SNR estimator probe and redirected the messages to message debug.)

I've tried the examples adapted to my input and M2M4 estimator spits out somewhat expected values (so I've been told by others in the company).

The code:

from sys import argv
from scipy import fromfile, absolute, array_split, mean, sqrt, log10, ceil, linspace, complex64
from matplotlib import pyplot as plt

def snr_est_m2m4(signal):
    M2 = mean(absolute(signal)**2)
    M4 = mean(absolute(signal)**4)
    snr_rat = sqrt(2*M2*M2 - M4) / (M2 - sqrt(2*M2*M2 - M4))
    return 10.0*log10(snr_rat)

def chunkify(lst, chunk_size):
    return [lst[i * chunk_size : (i+1) * chunk_size] for i in range(int(ceil(len(lst)/chunk_size)))]

def main():
    if len(argv) != 2:
        return
    bits = fromfile(argv[1], dtype=complex64)
    m2m4_axis = [absolute(snr_est_m2m4(s)) for s in array_split(bits, 100000)]

    f1 = plt.figure(1)
    s1 = f1.add_subplot(1,1,1)

    chunk_size = 100
    smooth_m2m4 = [mean(c) for c in chunkify(m2m4_axis, chunk_size)]
    s1.plot(linspace(0, 1000, len(smooth_m2m4)), smooth_m2m4, 'b-', label='m2m4')
    s1.grid(True)
    s1.set_xlabel('time')
    s1.set_ylabel('snr')
    s1.legend()
    print("Done!")
    plt.show()

if __name__ == "__main__":
    main()

And the output I receive:

I know the recording is valid since the GRC chain we currently have implemented gives out an image of the Earth.
I've tried the estimators on a different signal but yet again, the values didn't make sense; the SNR was higher when there was no signal.

What do the parameters alpha and Samples between SNR messages in GRC represent? Why don't the estimators in GRC work as expected? What is the suggested way of calculating SNR?

I'm at a loss here, please send help.


texts: Significance of an impulse in the frequency domain

texts:
I know that $X(f)$ gives the amplitude associated with the frequency component $f$ of a signal $x(t)$.

Now, a sinusoidal signal in time $x(t) = A \cos (2 \pi f_0 t)$, has a Fourier transform $X(f) = \frac{A}{2}[\delta(f-f_0) + \delta(f+f_0) ]$.

My question is that the Dirac Delta funcion tends to $\infty$ at $0$. Then, multiplying it by $A/2$ should also result $\infty$. If this was the case, then what is the significance of $X(f)$ for $x(t) = A \cos(2 \pi f_0 t)$.

More specifically what is the significance of the impulses (having infinite amplitudes) in the frequency domain?

Thank you. :)


texts: What is the BER for multilevel modulated signal

texts: Assume I want to transmit a low or high symbol which can be BPSK modulated. For ordinary BPSK, bit 0 can be mapped to 1+i0 and bit 1 mapped to -1+i0. In my question, each two bits are BPSK modulated and a level dependent so bit 01 is mapped to 2+0i and bits 00 is mapped to 1+0i. Also bits 10 is mapped to -1+0i and bits 11 is mapped to -2+0i. This is BPSK (not QPSK) signal since there are BPSK demodulator and level detector at the receiver in which each of them output one bit at symbol period.I want to know what demodulation scheme to use in order to detect the symbol level and hence calculating the BER.


texts: Bandpass general equation to difference equation

texts:
I know this is a very basic question and I am coming out from a quarter of DSP. I want to create a function in Java which can taken in two parameters, either (double centerFrequency, double bandWidth) or (double startFrequency, double cutoffFrequency).

Based on the information I've gathered, I want a non-ideal IIR filter because this is for a project I am doing with an accelerometer that needs to be updated real-time. Hopefully I am not going too much into programming, but the general form of my filter function would be like this:

double BPFilter(double x, double centerFrequency, double bandWidth) {
    double filteredY;
    static double[] y = {0, 0, 0};
    //y[0] represents y[n]
    //y[1] represents y[n-1]
    //y[2] represents y[n-2]

    y[0] = A*y[2] + B*y[1] + C*x;
    y[2] = y[1];
    y[1] = y[0];
    return y[0];
}

What would be the best way to calculate A, B, and C? Currently I know how to do continuous transfer functions, $H(s)$, so if that is translatable to A, B, and C I could do it that way. The sampling rate is at $100\textrm{ Hz}$.

TL;DR: Discretize a continuous transfer function by hand (no MATLAB)


texts: Understanding Arrangement of second order sections

texts:
When splitting higher order recursive filters into more convenient 2nd order sections (biquads), in matlab, we usually use zp2sos to convert. In document, there's a couple of lines , which I am quoting below,

[sos,g] = zp2sos(z,p,k,''order') specifies the order of the rows in sos, where 'order' is
'down', to order the sections so the first row of sos contains the poles closest to the unit circle
'up', to order the sections so the first row of sos contains the poles farthest from the unit circle (default)

Here i can't remember that I had read something relevant to arranging the biquads on the basis of poles location.
Now Can anyone explain me,
why should we arrange the biquads in a specific manner while cascading them ?
What that arrangement is to deal with the poles being closer or farthe from unit circle, as Matlab documents says?
What specifics need to be considered while ordering the biquads to get higher order filters.
It'll be so nice to illustrate by using some case study, for example, designing the A weighting filter (8th order IIR), and converting it ino biquads, and arranging them in a specific manner which is under discussion.
Thanks in advance.


texts: What is the exact meaning of unstable system in DSP?

texts:
In physical systems I understand what is the meaning of stability or unstability. An operational amplifier for example, if working in positive feedback will either saturate or start osscilating (i.e will not have any stable state). thats clear to me.

But I am unable to understand what exactly we mean when we say an IIR filter (or any other digital system), for example, can become very unstable.

What exactly inside the Digital Signal Proccessor  happens, what happens to the output physically?
What exactly we mean by unstable system in this context?


texts: Required number of measrments for signal recovery in a compressed sensing MMV problem?

texts:
For multiple measurement vector (MMV) problem
$Y=AX$
where $A$ is $m \times n$ sensing matrix and  $X$ is $n \times L$ matrix haveing K non zero rows.

What are the necessary conditions on the number of measrments $m$ and sparsity K to recover the signal $X$ from $Y$ using OMP?


texts: Incorrect Butterworth filter outcome using Rational Transfer Funcion

texts:
I've been playing with Second-Order Butterworth high-pass filter for a while.
I've calculated the coefficients using Matlab function and I wanted to use them to filter simple sine wave generated by Matlab.
Unfortunately the results are not correct/consistent, verified by Matlab filter(b,a, signal) function.
My Matlab code is as following (ommiting plotting part):

clear all;
fp = 10000; %sampling frequency
ts = 0:1/fp:0.05; %number of samples, length
squareF=100; %frequency of signal
squareWave = sin(2*pi*squareF*ts); %generate the sine
fc = 500;
[b,a] = butter(2 ,fc/(fp/2), 'high'); %numenator, denumenator
for n = 3:501
    filteredWave(n) = (b(1)+b(2)*squareWave(n-1)+b(3)*squareWave(n-2))/(a(1)+a(2)*squareWave(n-1)+a(3)*squareWave(n-2))*squareWave(n);
    %Rational Transfer Funcion
end
filteredCorrectly = filter(b, a, squareWave);

Filter plots are below.

Where is the error?


texts: Inverse Fourier Transform of a squared sinc function

texts:
I have here a squared sinc function, which is the Fourier Transform of some triangular pulse:
$$\mathrm H(f)= 2\mathrm A\mathrm T_\mathrm o \frac{\sin^2(2\pi f \mathrm T_\mathrm o)}{(2\pi f\mathrm T_\mathrm o)^2}$$
As an excercise, I would like to go back to the original time domain triangular pulse, using the inverse Fourier Transform.
I feel like I'm very close to achieving it, however, I stumbled upon this integral:
$$\mathrm h(t)=\frac {\mathrm A} {2\pi ^2 \mathrm T_\mathrm o}\int_{0}^{\infty}\frac{(1-\cos(4\pi f\mathrm T_\mathrm o))\cos(2\pi ft)}{f^2}df$$
I'm not sure if there is already an expression which defines the value of this integral, or if this integral is of a well-known type, with a well-defined value. The tables of integrals that I have right now are basic and don't cover integrals of such type.
Does someone know if there is a well-known expression for the value of integrals of this type?
Or, would you be so kind as to recommend a good table of integrals for Fourier Analysis?
Thank you very much!


texts: Block normalization of histogram of oriented gradient

texts:
I am a beginner in image processing. I have some questions regarding HOG.

https://www.learnopencv.com/histogram-of-oriented-gradients/

From this website, it mentioned that each element inside the HOG descriptor will be divided by the sum of all elements in a block. However, how to do block normalization if overlapping the blocks? I do not understand because some elements in the vector will be elements of more than one blocks.

Picture of block normalization

For example, if we have a histogram [v11 v12 v13 v14...v19] ,[v21 v22 v23 v24...v29], [v31 v32 v33 v34 v35...v39] and [v41 v42 v43 v44..v49] in the red block, we can normalize the first element of the first vector like this: v11/(v11+v12+...+v21+v22+v23+..+v31+v32+...+v41+v42...+v49)

However, how can I normalize the vectors in the red block and yellow block?

Thank you very much.


texts: Analytical Solution for the Convolution of Signal with a Box Filter

texts:
I have an exercise in which I am trying to filter an input signal $y(x) = \sin(x)$. Ideally, I would like to apply a box filter to this signal.

Previously, I successfully convolved the input signal $y(x)$ with a decaying response $h(x) = e^{-x}$.

I did so by the following the definition of convolution (e.g., integrating $\int_0^t\sin(x')e^{-(x-x')}\mathrm{d}x'$ and computing a damped sinusoidal signal.

My box filter is given by $\frac{1}{\Delta}$ for$|x-\xi| \leq \frac{\Delta}{2}$ and 0 elsewhere, where $\Delta$ is the filter width. I understand that a box filter is a local average, and I can implement this numerically, but I do not understand how to analytically integrate this as I did with the damped exponential 'filter'.

I tried to take the Fourier transform of $y(x)$ and $h(x)$ and multiplying them in Fourier space, but I could not figure out how to do so.

Thanks for any help.


texts:
I'm working on a research project where we would like to apply convolutional neural networks to an image representation of a signal. However, it seems that if I would use a spectrogram, I would end up loosing the phase information. As a result, my question is:

Does an equivalent transformation of a signal to a spectrogram image exist in which the phase information is part of the resulting image?

Thanks a lot and have a great day, Maxime

From https://en.wikipedia.org/wiki/Spectrogram:

  Limitations and resynthesis

  From the formula above, it appears that a spectrogram contains no
  information about the exact, or even approximate, phase of the signal
  that it represents. For this reason, it is not possible to reverse the
  process and generate a copy of the original signal from a spectrogram.


texts:
The Hilbert Transform of a 1D/real-valued vector signal returns the analytic signal, x, from a real data sequence, xr. The analytic signal x = xr + jxi has a real part, xr, which is the original data, and an imaginary part, xi, which contains the Hilbert transform.
hilbert uses a four-step algorithm:

1. Calculate the FFT of the input sequence, storing the result in a vector x.

2. Create a vector h whose elements h(i) have the values:

     1 for i = 1, (n/2)+1

     2 for i = 2, 3, ... , (n/2)

     0 for i = (n/2)+2, ... , n

3. Calculate the element-wise product of x and h.

4. Calculate the inverse FFT of the sequence obtained in step 3 and returns the first n elements of the result.

This algorithm was first introduced in [8]. A python implementation of this can be seen bellow:
from scipy import linalg, fft as sp_fft
import numpy as np

def hilbert(x, N=None, axis=-1):

        x = np.asarray(x)
        if np.iscomplexobj(x):
            raise ValueError("x must be real.")
        if N is None:
            N = x.shape[axis]
        if N <= 0:
            raise ValueError("N must be positive.")
        print(x.shape,N,axis)

        Xf = sp_fft.fft(x, N, axis=axis)
        print(Xf.shape)
        #plt.plot(Xf)
        #plt.show()
        h = np.zeros(N)
        #plt.plot(h)
        #plt.show()
        if N % 2 == 0:
            h[0] = h[N // 2] = 1
            h[1:N // 2] = 2
        else:
            h[0] = 1
            h[1:(N + 1) // 2] = 2
        print(h)

        if x.ndim > 1:
            ind = [np.newaxis] * x.ndim
            ind[axis] = slice(None)
            h = h[tuple(ind)]
        x = sp_fft.ifft(Xf * h, axis=axis)
        return x

This code was taken form Scipy implementation Scipy.signal.hilbert. I am looking to invert/reverse this process (inverse_hilbert), the best description i have found to do this is from Mathworks Inverse Hilbert Transform

However the real and complex arrays currently have me at a loss not sure how this feeds into this equation, or if this is the correct equation as wikipedia has a different equation to my understanding.

If we create a random complex signal and compute the hilbert Transform on it we get the following 2 arrays one real and one imagery. I am looking to reverse this transform any help would be appreciated.
import matplotlib.pyplot as plt

x = np.arange(0, 30, 0.1);
y = np.sin(0.05*x)+np.sin(6*x)+np.cos(3*x)
plt.plot(x,y)
plt.title('Complex Waveform')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

x_a = hilbert(y)

plt.plot(x,x_a.real, label='Hilbert Real', alpha=0.5, lw=2)
plt.plot(x,x_a.imag, label='Hilbert Imag', alpha=0.5, lw=2)
plt.grid()

Any comments here to help my under standing would be appreciated, really looking for a step by step way to reverse this transformation.
[8]: Marple, S. L. “Computing the Discrete-Time Analytic Signal via FFT.” IEEE® Transactions on Signal Processing. Vol. 47, 1999, pp. 2600–2603.


texts: When do we normalize data in audio processing?

texts:
Lets assume we have a simple audio processing flow with preprocessing, feature-extraction and classification.

Where and why do we normalize the data?

I know the purpose of normalizing the data in the preprocession-step. But as I implemented some feature-extracrion methods like spectral flux, I have read about normalization of the extracted features.

Do we need to normalize the extracted features, when we normalized the data in the preprocessing step and vice versa?


texts: Why does the solution for the voltage of this RLC circuit result in a maximum capacitor value?

texts:
I have a circuit with an inductor (L), capacitor (C), and two resistors (R1=R2):

This represents an equivalent electrical circuit to a plectrum or hammer (cap=spring, inductor=mass) striking a string divided into two segments (R1 & R2).

The system is energized by a current impulse at time 0. This DSP synthesis technique is explained more here: https://www.dsprelated.com/freebooks/pasp/String_Excitation.html

The problem I am having is that the solution for the voltage leads to a condition where there is a maximum capacitor value or it is unsolvable, and I am not sure why or what this represents.

The current across the parallel component (the "exciter") is described as follows:

$I(s) = I_C(s)+I_L(s)$

$I(s) = sCV_C(s) + \frac{1}{sL}V_L(s) + \frac{i_0}{s}$

$I(s) = \frac{CLs^2+1}{Ls} V_{exciter}(s) + \frac{i_0}{s}$

The current across each string is represented by:

$I(s) = V_R(s)/R$

Furthermore, if $R1=R2$, $V_R(s)$ can be approximated given that with the exception of the momentary impulse:

$V_{exciter}(s) + V_{R1}(s) + V_{R2}(s) = 0$

$V_{exciter}(s) + 2V_{R}(s) = 0$

$V_R(s) = \frac{-V_{exciter}(s)}{2}$

Therefore we can put these equations together as:

$\frac{-V_{exciter}(s)}{2R} = \frac{CLs^2+1}{Ls} V_{exciter}(s) + \frac{i_0}{s}$

$V_{exciter}(s) = \frac{-2i_0LR}{2CLRs^2+ Ls + 2R}$

To then convert this to the time domain, I need it in terms of: $\frac{1}{(s+a)(s+b)} = \frac{e^{-at}-e^{-bt}}{b-a}$

This is solved as:

$V_{exciter}(s) = \frac{-2i_0LR}{2CLR(s^2+s/(2CR)+1/(CL))}$

$V_{exciter}(s) = \frac{-i_0}{C}\frac{1}{(s+a)(s+b)}$

$a=\frac{-1}{4LR}\left(-\frac{L}{C}-\sqrt{\frac{L^{2}}{C^{2}}-\frac{16CLR^{2}}{C^2}}\right)$

$b=\frac{-1}{4LR}\left(-\frac{L}{C}+\sqrt{\frac{L^{2}}{C^{2}}-\frac{16LR^{2}}{C}}\right)$

This results in a final solution of:

$v_{exciter}(t) =\frac{-i_0}{C}\frac{e^{-at}-e^{-bt}}{b-a}$

The problem now is the square root term for $a$ and $b$. This creates a condition where in order to be solvable:

$L (L - 16 C R^2) >= 0$

So at a given $L$ and $R$ value, $C$ must be $< \frac{m}{16R^2}$ for this equation to be solvable or you otherwise end up with a square root of a negative which I am not sure how to deal with and get real results from.

Why does this happen and what does it mean? Is there some way to solve this equation differently so it will still work at any $C$ value?

It was suggested to me that above this capacitor threshold, the circuit may be going from exponential decay to sinusoidal.  But even still, should I not be able to work out an equation for the voltage that will still give me a value?


texts: Edge linking by adaptive mathematical morphology

texts:
I'm trying to implement the Adaptive Morphological Edge-Linking Algorithm (Frank Y. Shih, Image Processing and Pattern Recognition: Fundamentals and Techniques, section 5.5.2 OR Article link). The problem is I'm having trouble understanding how to implement the adaptive morphological dilation and how it should be applied (from the algorithm's description, I understand that I have to apply adaptive dilation at each end point, but wouldn't this affect the rest of the image too? I think I'm missing/misunderstanding something)

EDIT: Actually, I think I understand what it means by applying at each end point. Apparently, each point in the image has its own kernel or so. It does seem to say something about this in the operation's definition. I'll look more into it, hopefully I can find a way to implement it.


texts:
I am reading Proakis book "DSP using Matlab", 3rd edition.
I am reading chapter 8, section 8.3, p. 402, and I am confused regarding the equation of poles (roots of denominator of system function) eq 8.47.
How this equation has been derived from eq 8.46 especially highlighted line how the right most side of eq 8.47 $Ω_ce^{j(\pi/2N)(2k+N+1)}$ has been derived from left side $(-1)^{1/2N} jΩ_c$.


texts:
I encounter a situation that it is very beneficial to limit the signal strength below noise floor.

I know that the Shannon–Hartley theorem suggest that it is possible:

I searched over the internet but cannot find any direct pointer.

What keywords should I look for? Are there any compiled resources on the topic?

Thanks


texts: Extrapolate a 2D array using Fourier Transform

texts:
I need to extrapolate a given 2D array to a larger domain, keeping the spatial frequency. This is the original field:

(the data file in numpy npz format and a Jupyter notebook to plot it can be found here)

The horizontal size here is 14.864408108 (critical wavelength), and the vertical size is 14.864408108/sqrt(3)

I would like to extrapolate it to a square domain of size, let's say, ten times bigger than the original horizontal size, but fixing the number of grid points to a given number (256x256 for example).

Currently I use numpy.tile to do that:

import numpy as np
from scipy import interpolate

def enlarge(field,mult_x=5,mult_y=8):
    y = np.linspace(2*mult_y*Y.min(),2*mult_y*Y.max(),int(np.round(2*mult_y*field.shape[0]/mult_x)))
    field2 = np.concatenate((field, field[::-1, :]))
    field4 = np.concatenate((field2, field2[:,::-1]), axis=1)
    large = np.tile(field4,(mult_y, mult_x)) # 17 just because it gives a number close to 256*10
    less = large[::mult_x, ::mult_x] # skip every 10 rows/columns
    #print(less.shape)
    f = interpolate.interp1d(y, less,axis=0)
    y_new = np.linspace(2*mult_y*Y.min(),2*mult_y*Y.max(),less.shape[1])
    field = f(y_new)
    return field

And I get this figure:

Can I use Fourier transform to do that?


texts: Convolution sum vs auto-correlation problem

texts:
I had an exam one of these days and one of the questions was:

"Knowing that an auto-correlation estimator of a sinal x[n] could be defined by:

$$
R_{xx}[k] = \sum_{n=-\infty}^{+\infty}h[n].x[n+k] , k=0,1,...,n-1
$$
and that the non-periodic convolution:
$$
y[n]=x[n]*h[n] = \sum_{k=-\infty}^{+\infty}h[k].x[n-k]
$$
show how to calculate the autocorrelation with a convolution operation"

And here is my reasoning:

let's say that:
$$
h[n] = x[-n]
$$
Therefore, we would have some $y[n]$ for the following convolution:
$$
y[n] = x[n]*h[n]=x[n]*x[-n] = \sum_{k=-\infty}^{+\infty}x[k].x[k-n]
$$
So it would be nice to invert the signal $y[n]$, finding then a $y[-n]$. However, we will get the same result as before since convolution is a commutative operation:
$$
y[-n]=x[-n]*h[-n]=x[-n]*x[n]=x[n]*x[-n]=\sum_{k=-\infty}^{+\infty}x[-k].x[-k+n]
$$
Therefore, $R_{xx}[n]=y[n]=y[-n]$, since the convolution sum is symmetric with respect to zero.

Question:

Is my reasoning correct by formal mathematics?
The professor gave me a zero in this answer because he says that my assumption $R_{xx}[n]=y[n]=y[-n]$ is wrong because the correct calculation of the autocorrelation has a complex conjugate in one of the terms (so it might not be symmetric with respect to zero). I agree with it even though he has not placed this formalism in the autocorrelation definition. So, following the exam's definitions, is my reasoning correct?

Also, can I say that this will always be true whenever the signal is purely real?
And by consequence of the question above: can some of these be demonstrated for complex signals?


texts: Does the signal reconstruction error in compressed sensing using $l_1$ norm minimization depends on the amplitude of non-zero coefficients and their location ?


texts: Efficient Magnitude Comparison for Complex Numbers

texts:
Is there a more efficient algorithm (or what is the most efficient known algorithm) to select the larger of two complex numbers given as $I + jQ$ without having to compute the squared magnitude as

$$I^2+Q^2$$

In particular binary arithmetic solutions that do not require multipliers?  This would be for a binary arithmetic solution using AND, NAND, OR, NOR, XOR, XNOR, INV, and shifts and adds without simply replacing required multiplication steps with their shift and add equivalents (or what would be equivalent in processing steps). Also assume the number is represented with rectangular fixed point (bounded integers) coordinates (I, Q).

I am aware of magnitude estimators for complex numbers such as $max(I,Q) + min(I,Q)/2$, and more accurate variants at the expense of multiplying coefficients but they all have a finite error.

I have considered using the CORDIC rotator to rotate each to the real axis and then being able to compare real numbers. This solution also has finite error but I can choose the number of iterations in the CORDIC such that the error is less than $e$ for any $e$ that I choose within my available numeric precision. For that reason this solution would be acceptable.

Are there other solutions that would be more efficient than the CORDIC (which requires time via iterations to achieve accuracy)?

Determining Best Answer

There was incredible work done by the participants (and we still welcome competing options if anyone has other other ideas). I posted my proposed approach to ranking the algorithms and welcome debate on the ranking approach before I dive in:

Best Approach to Rank Complex Magnitude Comparision Problem


texts: Meaning of wavelet and scaling coefficients

texts:
What is the meaning of wavelet coefficients and scaling coefficients? E.g. for a sequence I obtained the following wavelet coeffients.

How am I supposed to interpret them? I used wavelets package in R. The sequence concers price time series of 21 days. I used the Haar wavelet.


texts: Detect Gradual Increase Before a Decrease in Noisy Data

texts:
I have data that'ss not too noisy and I am trying to detect a pattern where it gradually increases then decreases in a short period of time (20 ticks? it should be roughly similar per session but can vary slightly between different sessions). This pattern usually happens in pairs but not always.

So looking at the top plot of the data over time, the pattern I am after is at time 100 and 300. you can see the data is roughly stable, then gradually increases and decreases and stabilizes, then after a while again increases, decreases and stabilizes.
As you can see there are also instances where it suddenly changes (550 & 800), or it increases then decreases over a longer period of time (700-800) but I need to detect my short pattern only.
What I have done sofar:
I have taken the variance of the last 10 values which is the middle plot. It gives me a relatively good idea when the data starts to deviate but, variance ignores the sign so I cannot tell if the data is increasing or decreasing and it doesn't differentiate between sudden and gradual changes.
So I decided to take the squared difference between each value and the previous value and keep the sign. I took the square to exaggerate the difference to make it more distinguishable from smaller variations and kept the sign to help identify when it increases then decreases.
I can't say
if diff[i] > diff[i-1] for 5 ticks, and then diff[i] < diff[i-1] for 5 ticks because diff is noisy & changes its sign to positive midway at 130 even though overall it was decreasing.
I feel like there would be a better way than trying to filter the noise and then applying that if statement above. But I'm drawing blanks at the moment.
Whats the best way to programatically detect a short but gradual increase then decrease in noisy numbers?
I've added a csv with the data at this pastebin: csv_pastebin


texts: Determining the pattern direction of a bitmap image

texts: How can I obtain average direction of the pattern shown in the figure below. It is the direction of the red arrow relative to the yellow (horizontal) line. Any ideas for an approach? I couldn't figure out a way to approach. This is a spatio-temporal image created from a video. Thank you.


texts: chirp z-transform for different output sizes

texts:
I am attempting to use the chirp z-transform for an application that requires arbitrary FFT output sizes less than or equal to the length of the input signal.

However, I've encountered an issue where the total output energy of the CZT for output lengths different from the input length varies considerably from the input energy.  Initially I tested this by windowing (using flattop) the input data prior to performing the transform.  The returned amplitude values where not consistent depending upon the frequency of the input tone.

I removed the window and repeated the experiment to find the total energy values were correct if the input tone was aligned on the bin center of the first FFT in the CZT transform, but wildly incorrect otherwise.

My best guess is the scalloping loss in the first FFT of the algorithm is somehow causing energy to not be preserved elsewhere.

It would seem to indicate that the CZT algorithm violates Parseval's theorem for output sizes different from the input size.

Does anyone have any suggestions as to what is wrong here?


texts:
I have a signal $x(t)$ for which I want to find the Nyquist frequency :
$$ x(t) = \frac{\sin{\pi t/2}}{\pi t/2} \ast \sum^\infty_{n=-\infty}\delta(t-10n)$$

I am trying to solve this in the time domain like so :
$$ \begin{split} x(t) &= \frac{\sin{\pi t/2}}{\pi t/2} \ast \sum^\infty_{n=-\infty}\delta(t-10n)\\
 &=\frac{\sin{\pi t/2}}{\pi t/2} \ast[...\delta(t+20)+\delta(t+10)+\delta(t)+\delta(t-10)+\delta(t-20)+...]\\
&=\frac{\sin{\pi (t+20)/2}}{\pi t/2}
+\frac{\sin{\pi (t+10)/2}}{\pi t/2}
+\frac{\sin{\pi t/2}}{\pi t/2}
+\frac{\sin{\pi (t-10)/2}}{\pi t/2}
+\frac{\sin{\pi (t-20)/2}}{\pi t/2}
+...\end{split}$$

Now shifting will not change the frequency so :

$$\begin{split} \omega_m &= \pi/2\\
2\pi f_m &= \pi/2\\
f_m &= 1/4\\
F_{\text{Nyquist}} &= 2\times 1/4 = 0.5\end{split}$$

But when I solved this quesiton in the frequency domain I obtained $0.4$ which is the correct answer.

What is the mistake in this method?


texts:
I've got a vector of $100 000$ numbers. All numbers are equal ($7000$ for example).

If I perform FFT over this vector, what will I get?
From my understanding, I should receive a fixed DC line. Is this true?
Is there any website (i've searched) which can simulate that?

Because I did tried to run FFT with Python np.fft.rfft and np.fft.rfftfreq and didn't get it...

What am I missing here?

(My question is derived from a broader perspective, so bear with me...)


texts: scipy.signal.decimate vs array operations

texts:
The following program generates a signal with a bit of noise according to some bits.  I then try to decode the signal using convolve and decimate.  2 problems I notice:

The whole this gets thrown out of wack depending on the 'type' parameters that I pass to convolve, is there a better way of lining up the signal?
The decimate function seems to mess things up - everything works better if just use array indices to access the data I want.  Am I doing something wrong with decimate?

Here is a picture of the output:

import matplotlib.pyplot as plt
from math import sin, pi
from numpy import ones, zeros
from numpy.random import random_sample as sample
from scipy.signal import decimate, convolve

N = 100
amplitudes = [1, 0, 0, 1, 0, 1, 1, 0, 1, 0]
pulse = ones(N)

#generate the transmit signal
tx_pulse = pulse + [0.2*sin(n) for n in range(1,N+1)]
signal = zeros(N * len(amplitudes))
for i in range(0, len(amplitudes)):
    for j in range(0, N):
        signal[i*N+j] = amplitudes[i] * tx_pulse[j]
signal = signal + 0.01*sample(len(signal)) #add noise

#now try to get the data back...
signal2 = convolve(signal, pulse/N, mode='valid')
signal3 = signal2 > 0.5 #bit decision

#dec = decimate(signal3, N, ftype='fir') #THIS MESSES UP THE THE FIRST SYMBOL
dec = signal3[0:N*len(amplitudes):N]

plot4 = zeros(N*len(amplitudes))
for i in range(0, len(amplitudes)):
    plot4[i*N] = dec[i]

print dec

plt.plot(signal)
plt.plot(signal2)
plt.plot(signal3)
plt.plot(plot4, 'ro')
plt.show()