texts: Differences between matlab's 'impulse' and 'lsim' command; create impulse response

texts:
I've got following problem:

assume you've this continous system 1st order:

$$ H(s) = \frac{1}{1+s} $$

s = tf('s');
sys = 1/(1+s);

I'm going to discretize this system via impulse invariant transformation and plot the results.

sysD = c2d(sys,1e-3,'impulse');
impulse(sys);
hold all;
impulse(sysD);

So this seems to work.

If I try to create the impulse response for 7 seconds manually I use:

imp = zeros(1,7*1e3);
imp(1) = 1;
t = 0:1e-3:7-1e-3;
y = lsim(sys,imp,t);
plot(t,y)

Unfortunately this is not the response I expected. The gain seems to be scaled by the samplerate 1e3.

So the question is: What does impulse do since this response should be correct. lsim uses the filter command and simply passes the nominator and denominator to filter I think. Is the dirac impulse (or kronecker delta) wrong? Must it be scaled by the samplerate?

I discovered the whole thing, because I wanted to design a very simple digital filter. I knew the responses of time continuous systems, so I thought it would be easy to discretize the transfer function and use the coefficients of the digital transfer function that can be used by filter or in a dsp application.

What is wrong?

Thank you very much!


texts: Position vector of Lukas-Kanade optical flow algorithm


texts: obtaning discrete-time system impulse response $h[n]$ from step response $s[n]$

texts:
The discrete-time system of figure below a step response: $$s[n] = \left[ 5 - 4 \left( \frac{4}{5} \right)^n  \right] u[n]$$

Find the impulse response of the system and verify it from the system
difference equation.
sr i want ans of this question.with explanation


texts: Difference equation and Impulse response of an IIR filter - deconvolution

texts:
Given the input to the IIR filter

$$
x[n]=\delta [n]+0.25\delta[n-2205], \quad \text{where}\quad \delta[n]= \begin{cases}1, &n=0\\0, & n\neq 0\end{cases}
$$

What I want to do is to design an IIR-filter that cancels the "echo", such that I'm only left with the $\delta[n]$. To do so I of course make an IIR-filter, where you take the $\mathcal Z$-transform of the impulse response, and cascade it with the inverted

$\mathcal Z$-transform of the input:

$$
X(z)=1+0.25z^{-2205}
$$

Transfer function of the IIR filter:

$$
H(z)=\frac{1}{1+0.25z^{-2205}}
$$

Then our output naturally is:

\begin{align}
Y(z)&=\left(1+0.25z^{-2205}\right)\cdot\frac{1}{1+0.25z^{-2205}}\\
& = 1\\
\Rightarrow y[n]&=\delta [n]
\end{align}

So the transfer function of my filter is $H(z)$. What I want to do is to find the difference equation for my filter, and possibly find the mathematical form of the impulse response. Since the order is so high I'm having difficulties.

I'm pretty sure the difference equation will not be infinite, while the impulse response will... why??


texts: Discrete Fourier transform - Norms of complex input signals and their transformation

texts:
Given a signal $\mathbf{z} \in \mathbb{C}^n$ and its Discrete Fourier  transform $\hat{\mathbf{z} }$, does $||\mathbf{z}|| = ||\hat{\mathbf{z} }||$ hold?

The question is given to me like this with no additional details. Information about what kind of norm is also not given. Does anyone have an idea what the question might be looking for?


texts:
I am attempting to compute the sum of dirac functions in the frequency domain via the following formulas:

$$ ReLu(f(x)) = f(x) * \sum_{i} \delta(x-x_i) $$
$$ x_i : f(x_i) > 0 $$
$$ F(ReLu(f(x)))(k) = (\sum_{i} e^{2\pi jk x_i})F(f(x)) $$

From what I understand I need to find the values where x is larger than zero, then sum $e^{2\pi jk x_i}$ over these values. I made myself a trivial example to determine if I am doing this correctly and there is no difference to the image after the convolution by the sum of the exponentials. What am I interpreting wrong here?

Note this idea is from the paper: http://cs231n.stanford.edu/reports/2015/pdfs/tema8_final.pdf

import numpy as np
import matplotlib.pyplot as plt

# Make image array with a single negative value
img = np.array([[-1.0,2.3],[5,7.8]])

# Dirac values
N = 2
x = np.arange(0, N, 1)
y = np.arange(0, N, 1)
xm, ym = np.meshgrid(x, y)
x0 = xm
y0 = ym
shiftMat = np.exp(1j * ((2.0 * np.pi)) * (xm*x0 + ym*y0))

# Set equal to shift mat
# In this trival example I know that [0,0] is only negative position
# So set to 0 and compute sum of all positions in which f(x) > 0 as detailed in paper
freqRelu = shiftMat
freqRelu[0,0] = 0
freqRelu = np.sum(freqRelu)

# Fourier Convolution and IFFT
imgFFT = np.fft.fft2(img)
freqR = np.multiply(imgFFT,freqRelu)
# freqR = np.outer(freqRelu,imgFFT)
reluedFreq = np.real(np.fft.ifft2(freqR))

# Spatial Relu For Comparision
reluedImg = img
reluedImg[0,0] = 0

plt.subplot(121)
plt.imshow(reluedImg)
plt.subplot(122)
plt.imshow(reluedFreq)
plt.show()

# print(np.allclose(reluedFreq,reluedImg))
print(reluedFreq)
print(reluedImg)


texts: How can I mimic this ADC for a given resolution, sampling rate and input with MATLAB or Python?

texts:
I'm quite new to the area so I need help to progress at some point.

Basically I have a signed 16 bit ADC which means it maps -10V to +10V range to 16 bit. So the max voltage +10V is 2^15-1 = 32767; zero is 0; and -10V is -2^15 = -32768.

I want to create/plot a sinusoid in MATLAB or Python where I can mimic exactly the ADC which samples a 1V 100Hz sinusoid input for 1 seconds which is sampled with 512Hz sampling rate. At the end I want to obtain the same FFT plot where I would obtain after a real ADC sampling.


texts:
One of the properties of convolution is that the centroids of the functions being convolved should add. See equation (18) of Wolfram Alpha's section on convolution Convolution Properties.

For example, if we have a Gaussian centered at t1=300, and we convolute it another Gaussian centered at t2=200, the resulting Gaussian peak after convolution is centered at 500. MATLAB works fine for that.

I was trying to convolute Gaussian functions, one centered at 0 (orange), and other one at 500 (blue), but the resulting function is not centered at 500 in MATLAB but appears much later (yellow). How can one circumvent this problem in MATLAB? I want to keep one Gaussian at 0. Thanks.

'''W=[0:1:2500]'; % bin numbers

g1= normpdf(W, 500, 100); % First Gaussian with unit area (x, mu, sigma)

g2= normpdf(W, 0, 100); % Second Gaussian with unit area (x, mu, sigma)

Convg1g2=conv(g1,g2)

plot ( g1)
hold on
plot( g2)
hold on
plot(Convg1g2)'''


texts: how to plot magnitude vs frequency of ECG signal taken from database in MATLAB?

texts:
load 100m
x= val(1,:);
fs = 1000;
Fn=fs/2;
Ts = 1/fs;
t = (0:length(x)-1)/fs;

subplot(3,1,1)
plot(t,x)
title('input signal')

y = fft(x);
fr=freqz(y);
f = (0:length(y)-1)*5/length(y);
subplot(3,1,3)
plot(f,abs(y))
title('Magnitude')


texts:
I have magnetometer time-series measurements in mV. I am trying to understand what is theory behind software that the company has provided.

Let $x$ be your time series vector, which is an array of double, and $W$ your
window length. $x$ is sliced into $L$ windows, each window is detrended and
(if not rectangular window function is chosen) a Hanning window is
applied. Then the FFTW is used to obtain an array of complex numbers $X$.
From $W$ real input values we get $P = W/2$ complex numbers in $X$. $X$ is calibrated with

$$
    X'= \frac{X}{\sqrt{f_s*P}}
$$
Transfer function is given by
$$
   H(\omega) =0.8V/nT*F_m(\omega)*F_t(\omega)*V
$$
where
$$
F_m(\omega)=\frac{j\omega}{j\omega+2\pi 4}\\
F_t(\omega)=\frac{2\pi 1591}{2 \pi 1591 + j\omega}
$$

If a transfer function $H$ is given we obtain

$$
    X'' = X'/H
$$

Finally if you use a windowing function the result has to be multiplied with the inverse of the integral of that window.

My question is why do I need the inverse of the integral of that window?

Are my units going to be $nT$?

I am giving part of the calibration file

 FREQUENCY    MAGNITUDE      PHASE
    Hz                     deg
Chopper Off
+4.0000E-01  +1.0714E-02  +1.4808E+02
+5.5580E-01  +1.3325E-02  +1.3911E+02
+7.7230E-01  +1.5775E-02  +1.2956E+02
+1.0731E+00  +1.7732E-02  +1.2025E+02
+1.4910E+00  +1.9094E-02  +1.1174E+02
+2.0718E+00  +1.9931E-02  +1.0442E+02
+2.8787E+00  +2.0392E-02  +9.8155E+01
+4.0000E+00  +2.0603E-02  +9.2520E+01
+5.5580E+00  +2.0611E-02  +8.7089E+01
+7.7230E+00  +2.0408E-02  +8.1384E+01
+1.0731E+01  +1.9949E-02  +7.4971E+01
+1.4910E+01  +1.9006E-02  +6.7315E+01
+2.0718E+01  +1.7761E-02  +5.8714E+01
+2.8788E+01  +1.5629E-02  +4.9042E+01
+4.0000E+01  +1.3116E-02  +3.9247E+01
+5.5580E+01  +1.0350E-02  +3.0040E+01
+7.7230E+01  +8.0188E-03  +2.2525E+01
+1.0731E+02  +5.9868E-03  +1.6408E+01

It seems that people do not understand my questions.Manuel can dowloaded from here

https://geo-metronix.de/ocl/public.php?service=files&t=9b54b82730e4b8b869e083574923d99a

You can go to tsmp and download tsmp.pdf.Calibration is on page 40.


texts: Conceptualising the continuous time unit impulse function as derivative of unit step

texts:
This a very newbie question.
I just watched Lecture 3 of Oppenheim's Signals course and he defines  here the continuous time function as the derivative of the unit step function like so:
$$ \delta_\Delta(t)=\frac {du_\Delta(t)}{dt}$$
and that $ \delta(t)= \delta_\Delta(t) $ as $\Delta \to 0$
He claims that the derivative is equal to 1 no matter the value of $\Delta$, because that derivative can be interpreted as the area of a rectangle with sides $\Delta$ and $\frac 1 \Delta$
I can't conceptualise this the way the function $u_\Delta(t)$ is drawn at all. If the function is linear, that is, $y = mx + b$ passes through the origin, meaning $b=0$, and we can see that it has the point $(\Delta, 1)$ we can easily tell that $m= \frac 1 \Delta$ and that should be the derivative.
Can someone explain to me the error in my line of thought? Why is the derivative the area of a square?


texts:
I have one doubt according to modulation of a real signal.

there are so many examples where the real signal was modulated with frequency bigger then 0.

Why the signal could not be putted in f=0?

is there connection with a sending power?
( Antenna send bandpass signal $\rightarrow$ signal is matched $\rightarrow$antenna send low P $\rightarrow$power will be reflected, Tx get P again)


texts: Algorithm for finding length and end points of a curly human hair

texts:
I'm trying to do some analysis of a curly human hair on an iPhone app. I've been using the openCV library.

Particularly I want to be able to estimate the length of the hair, and locate its end points (a measure of curliness is stretched end-to-end length vs. relaxed end-to-end length).

My approach so far has been to use a canny edge detector to find the hair. Then to dilate the result. Finally to skeletonise the image and consider the ends to be those pixels that only have one neighbour. Since the resulting image would hopefully be a one pixel thick line, I could just count the pixels to estimate the length (and maybe adjust for diagonal movements).

The trouble is when the hair is so curly it crosses itself. This makes the skeletonisation add in extra branches and rings that mess up the results.

Is there a better approach for very curly hairs?

Some examples:


texts:
An ISI channel of length 'L+1' $\mathbf{h}=[h_0 \space h_1 \space … \space h_L]$can always be modeled as a finite state machine(FSM) of memory 'L'. For a binary input $x\in\{+1,-1\}$,a state can be defined as $\mathbf{s}=[x_{1}\space  x_{2} \space … \space x_{L}]$.

I recently came across a paper (cf. this article):

Myint, L.M.M.; Supnithi, P.; Tantaswadi, P., "An Inter-Track Interference Mitigation Technique Using Partial ITI Estimation in Patterned Media Storage," Magnetics, IEEE Transactions on , vol.45, no.10, pp.3691,3694, Oct. 2009

in which $\mathbf{h}=[h_{k-1} \space h_k \space h_{k+1}]=[\beta \space b \space \beta]$, where the authors define gain=$b$ and ISI-coefficient=$\beta$ (ref:eq 3,4,5) . the noiseless channel output is defined as

\begin{equation}
y_{k} = bx_k+\beta(x_{k-1}+x_{k+1})\\
r_{k} = y_k + n_k
\end{equation}

and create a 4-state trellis using $[x_{k-1} \space x_{k+1}]$.

The branch metric is computed as $\lambda_k=(r_k-y_k)^2$.

This is approach is different from the conventional state definition $[x_{k-1} \space x_{k}]$. Since this paper is published in IEEE transactions on magnetics, I am hesitant to say that there is a mistake here but on the other hand this definition doesn't sound right to me.

I will be very grateful if someone can comment on this.

$\color{red}{New \space Edit:}$
While searching for the answer I came across another paper
(cf. this article)

Rubsamen, M.; Gene, J.M.; Winzer, P.J.; Essiambre, R., "ISI Mitigation Capability of MLSE Direct-Detection Receivers," Photonics Technology Letters, IEEE , vol.20, no.8, pp.656,658, April15, 2008
doi: 10.1109/LPT.2008.919597

That takes a similar approach.


texts:
Given a signal $$x(t)=16\cos(20\pi t+\frac\pi 4)+6\cos(30\pi t+\frac\pi 6)+4\cos(40\pi t+\frac\pi 3)$$how can I calculate the power contained in a frequency interval, say 12Hz to 22Hz.
The total power is $\frac{16^2}{2}+\frac{6^2}{2}+\frac{4^2}{2}=154$


texts: Difference Equations, why M<=N for causality?

texts:
in my notes for DSP they have the difference equation general form as:

y(k) + a1*y(k-1) + a2*y(k-2) + ... + an*y(k-N) = b0*x(k) + b1*x(k-1) + ... + bm*x(k-M)

with the claim that for the output to not depend on future values of the input then M<=N.

Why is this the case?

In fact, if I'm not missing something here I think I can provide a counterexample:

y(k) = x(k-1)

which has N = 0 and M = 1.
But this is just a unit delay, the output is just the previous value of the input.

What's going on here?


texts: Which book of Proakis has dealt with continuous time signals and systems?

texts:
I have heard of three books of Proakis namely Communication System Engineering, Digital Signal Processing and Digital Communications. I love the way he writes and would like to read about continuous/analog signals and systems by the same writer. Some of it has been dealt with in his book 'Communication System Engineering' namely Fourier Series and Fourier Transform. I want to read other continuous time topics like Laplace transform by the same author.

  Is there any book by Proakis available which fully deals with
  continuous time signals and systems?


texts: Anti-aliased boolean masks for integration of CCD images

texts:
(This question was already asked on Stack Overflow under the tags Matlab/Octave/Signal-processing/Image-processing, but I am reposting it here following the advice of one user that this site may be a more appropriate forum for such questions.)

Is it common/useful to generate an anti-aliased boolean mask for CCD image integration purposes ?

Let's imagine that I want to integrate the light that strikes a given area of a CCD sensor. The boundaries of this area correspond to physical coordinates that have no obligation to coincide exactly with the pixels of my CCD sensor.

In Matlab/Octave/(or even Scientific Python), the common algorithm used for integrating such an area is to define a boolean mask using a logical operation on an array, such as:

mask = R < radius;

(where 'R', and therefore 'mask', are 2D arrays and 'radius' is a float). Such a mask has values of 0 or 1.

I can then integrate the pixels that are comprised within the boundaries of my area by summing the masked image:

integrated_signal = sum(sum(mask.*image));

(where 'image' is the output of the CCD sensor, a 2D array with U16 values for instance).

However, mathematically speaking, nothing would prevent me from defining a continuous-valued mask, that is, a mask whose pixels could take any value between 0 and 1. Physically, this means that I could integrate portions of pixels. I could even compute such a mask using an anti-aliasing algorithm in order to create better approximations of non-square-friendly shapes such as circles (for integrating a signal with revolution symmetry, let's say).

My question to digital imaging geeks is the following: Is it standard practice to define such anti-aliased mask for integration purposes of imaging sensors ? Is it useful (as I suspect), or does it produce physically irrelevant results ?


texts: Calculate the Standard Deviation of Fundamental Frequency (MFCC)

texts:
I'm implementing a gunshot detector following the article "Algorithm for Gunshot Detection Using
Mel-Frequency Cepstrum Coefficients
(MFCC)" (paywall).

In the article, the authors uses 22 features based on MFCC (coeficients normalized between 0.1 and 0.9). One of the features is "Standard deviation of fundamental frequency". I already spent all my day trying to figure out the proper way to do this calculation.

The expected value for a gunshot is near 0.04.


texts: Is it possible to calculate/estimate the storage space a certain duration of .mp4 video files will take up? I have the resolution (eg. 4000 X 3000), the format (H.264 or H.265), the audio format (AAC) and the sample rate (41 000 Hz).


texts: How long does it take for snd_pcm_writei() (ALSA) to send data?

texts:
Using a pretty typical call to Linux's ALSA, snd_pcm_writei(pcm_handle, pcm_buffer, pcm_buffer_size) how do you know how long it should take to send this data? What is the formula?

The sound card should send frames (a frame is two samples for stereo, in ALSA terminology) at a certain rate, however, what is this rate? Sample rate? I don't understand how the sample rate relates to this.

If I have a sample rate of 44100, does it mean it sends 44100 samples per second? Given that my buffer in my program contains 8820 samples, I should see it drain this buffer 5 times every second (44100 / 8820 = 5). However, in reality I see it drain the buffer something like 10 times per second.

I could provide code, but I think it only would make the question less clear.

I've read https://stackoverflow.com/questions/24040672/the-meaning-of-period-in-alsa and www.linuxjournal.com/article/6735?page=0,1#N0x19ab2890.0x19ba78d8 as well as some other guides. However, it is still unclear to me how long it actually takes for the sound card to send data that it has received from snd_pcm_writei().


texts:
A system given by $\frac{s-1}{(s+1)(s-2)}$ has to be inverse transformed so that it is anticausal and nonstable. The answer given is $h(t)=-\frac{1}{3}(2e^{-t}+e^{2t})u(-t)$

Why the minus sign at the beginning?


texts: Using neural networks for rotation invariant OCR

texts: I have video with text which can be rotated in any direction. If I trained a neural network on images of each letter at all possible degrees of rotation, would it be possible to classify the characters accurately? Alternatively would another technique like HAAR, or KNN be better?


texts:
I am a PhD student with little engineering experience, currently attempting to recreate the signal processing techniques used in this paper: Monitoring Vital Signs and Postures During Sleep Using WiFi Signals. As a TL;DR: changes in channel state information (amplitude, in this instance) can be attributed to movements in otherwise still environments, most notably breathing and heart rate.

In recreating the signal processing performed in this paper, I have found it challenging to follow how the authors produced results as seen in Figure 8a and 8b, which should show isolated breathing frequencies.

At this phase, the signal (across 30 subcarriers) has been filtered using:

Hampel filter (for significant outlier removal)
Running mean filter (for smoothing)
Butterworth bandpass filter (to isolate potential breathing frequencies)

Once this has been done, a heatmap of amplitude across the subcarrier spectrum can be produced which should show clear banding, as seen in Figure 8a.
The Power Spectral Density for each subcarrier can be plotted, with the peak of the mean indicating the strongest observed frequency (Figure 8b). However, the lines plotted in this figure are curved, and I have been unable to find another example of curved PSD. It seems as though there is less data than implied in this graph, as this figure represents 5s of data (which was captured at 20Hz), which means the PSD could only be plotted for a 100n FFT. In my understanding, this should mean the bins are 0.2Hz wide? It appears as though they are able to use these curves across the peaks to infer additional information.

Using MATLAB's spectrogram, I have been able to plot the observed PSD in a similar signal (captured at 100Hz), and the peaks on this graph are not curved. My understanding of FFTs is that the peaks for each bin can be observed, but an asymmetrical curve as seen in Figure 8b cannot be drawn as we have no further information other than the bin.

Is my understanding of FFTs or PSD incorrect, or have I missed an obvious step in producing the results as seen?


texts: White noise generation for different bandwidth signals in MATLAB

texts:
I'm working on simulating an uplink scenario in satellite communication. The goal is to send a preamble + data to the satellite from ground station. Data symbols have a bandwidth of $15\ \rm kHz$ and the preamble has a bandwidth of $180\ \rm kHz$.
I'm simulating the transmitter, channel and receiver in MATLAB. I've a question regarding how to introduce white noise.
I read that variance of white noise depends on the bandwidth of the signal.
$$SNR = \frac{P}{NB}$$
Where $P$ is the power of signal, $N$ is PSD of white noise, $B$ is the bandwidth of the signal. In above case, the signal contains 2 parts, preamble + Data. Preamble has $180\ \rm kHz$ of bandwidth, whereas data has only $15\ \rm kHz$.
How to generate the white noise in this scenario? What should be the variance of the noise?


texts: Is subtraction a common method for efficient filtering?

texts:
I've recently been doing some DSP programming with regard to filter banks. One topic that I've come across was that of the dyadic filter bank. This filter bank essentially breaks a signal into sub-band signals corresponding to some octave.

Looking at some diagrams, I notice that many of them illustrate implementing both a highpass filter and low pass filter. This will essentially generate two sub-band signals. Conceptually, this makes sense.

However, it seems rather inefficient to implement two separate filters. Would it not be more efficient to implement a single filter, and then subtract the output from the original signal to determine the other signal?

In other words, implement just the low pass filter to obtain the lower subband signal. Then simply subtract that lower subband signal from the inputted signal. This should yield the higher band signal.

So my question:

Is this technique common?
Furthermore, would it be vulnerable to phase distortion in the sense that your low pass filter may have phase distorted or lagged the lower band signal?


texts:
I have a basic implementation knowledge of the 2D 8x8 DCT used in image & video compression. Whilst reading about Principle Component Analysis, I can see a lot of similarity, albeit PCA is clearly more generic. When I've read about DCT previously it was always presented in relation to DFT. So my question is how can the DCT be derived from a PCA perspective? (even a hand-wavey explanation is sufficient)

Many thanks


texts: Math for converting downward compressor to upward expander

texts: I’m looking to make an upward compressor.  Could it be as simple as making the compressor ratio a negative number OR making the numerator less than 1?


texts: How do I construct input to neural network from audio signals?

texts:
Input: Microphone recordings of digits from 0 to 9 from different speakers.

Output: The digit from 0 to 9.

I am doing this for fun. So first I will train my neural network using some samples and then use it to classify digits.

Problem is every person takes different amount of time to say the digit also every person takes different amount of time to say different digits. This gives me sound signal of different time duration. But input to the neural network  must be of some fixed size. How can I make all my input signals of same length so that I can give input it to neural network? Is there any standard technique? Fourier transform?


texts: Number of datapoints after decimation in discrete wavelet wavelet analysis

texts:
I am following an example provided by MATLAB on 1-D wavelet decomposition (DWT). The number of datapoints is 10000. If we do a level 3 decomposition following a similar example provided by MATLAB using https://www.mathworks.com/help/wavelet/ref/wavedec.html

The length of the coefficients "l" after wavelet decomposition is

[127 127 252 501 1000], corresponding to the length of approximation coefficients at level 3, number of detail coefficients at level 3, number of details coefficients at level 2 and number of detail coefficients at level 1. Why are these numbers not exactly divided by 2?

For example, 1000/2 = 500, it is 501. Similarly at next level 500/2 should be 250 but it is 252 and at the third level 250/2 should be 125 but it is 127.

Secondly, every textbook says that DWT needs 2^n data points. Even if we choose a number such as 8192, the number of decimated coefficients are not exactly and successively divided by 2. What is the reason behind this?

Thanks.


texts: Find beginning of signal for BPSK Demodulation

texts:
I want to store a small amount of data in a sound file recorded by a microphone but I know nearly nothing about signal processing. I managed to demodulate my data encoded using binary phase shift keying with this schema:

(From here: http://courses.cs.washington.edu/courses/cse561/10sp/project_files/FinalReport_FCI_YH.pdf)

To find the beginning of the signal I tried searching for the time where the strength of my carrier frequency in the signal exceeds a certain threshold but this is way too unreliable.

What would be a good way to solve this?


texts:
I am trying to understand and visualize the concept of a covariance matrix. Suppose I have a matrix:

$A = \begin{pmatrix} 2 & 3 & 4\\ 5 & 5 & 6 \end{pmatrix}$

how do I calculate its covariance matrix, and what is its practical meaning?
(All I was able to understand by now is that on the diagonal of the covariance matrix, variances for particular variables are placed, and on the upper and lower fields are the correlations between those variables.)


texts:
I've recently asked a question on stack overflow asking for different topics one can visit to solve the problem of wave re-construction. One of the recommended solutions was to implement the Hilbert-Huang transform, as it can help in understanding non-linear/natural data.

Here are some samples and details surrounding their capture to help me state my question.

I am collecting pressure readings from a pressure sensor at 30Hz. When I put pressure on the sensor, the Y axis above expands. When I release the pressure, it slowly approaches towards the idle state. In this case, the idle state has a base pressure reading of roughly -30.

Samples ranging from 0 to 200, and even those ranging from 400 to 750, have a normal pattern of pressure readings. They are very natural, in that they are not always the same when increasing and decreasing. They are also very natural in that their pressure idles around -30.

The same cannot be said about the samples ranging from roughly 270 to 370. The reason that the samples looks significantly different than the remainder of the sampled data, is because I bend the sensor as I am applying and releasing pressure. Instead of idling around -30, my sample readings drop down to as low as -140. Note: I have highlighted the area of impact in red in the below image.

Ideally, what I'd like to play around and perform is re-construct that wave to be in line alike the rest of the data. The green line in the below graph draws out what i expect to see as a final result of re-constructing the wave. ( It's just a rough sketch to explain the idea )

With that in mind, onto what my question is about.

Below you will find the generated HHT results of my sample data, with MATLAB.

Interestingly enough, as soon as IMF of the 2nd and 3rd round come about, the "abnormal" noise is very evident between samples of 300 to 350ish. Ideally, what I'd like to do is minimize this tone to be in line with the rest of the data.

Questions:

While my noob HHT understanding is poor, would it be logical to set the stop criteria ( in this case ) to IMF of C2 or C3, and use the results to somehow subtract out the IMF from the original wave to re-construct the wave as i desire?

Is there maybe a better way I can expand my utilization of HHT to solve my problem? ( Any topics or brainstorming ideas would help )

UPDATE after Cye's reccoemndation
I managed to filter out CF2 and CF3, however, the results were not as I expected. ( Image below ) The graph only became smoother, but that large dip was still present.


texts: Checking if u[n] and v[n] are orthogonal

texts:
Let $x[n]$ and $y[n]$ be two orthogonal baseband signals, with bandwidth B. Check if the following signals are orthogonal:

$u[n] = x[n]*\cos(2*\pi*f_c*n)$

$v[n] = y[n]*\cos(2*\pi*f_c*n)$

I'm not sure what the problem asks to be honest.
1) What exactly does he mean by bandwidth?? I only know the definition of a bandwidth of a sinc function, but how can I think of it generally?

2) I assumed the case for a rectangular function and reduced the problem to:

$\sum_{n=-\frac{1}{B}}^{n=\frac{1}{B}}\frac{1}{2}\cdot x[n]\cdot y[n]\cdot \cos(4\cdot \pi\cdot f_c\cdot n)$

But I don't know where to go from there. It's also weird for me to evaluate since $\frac{1}{B}$ is likely not an integer.

EDIT: also, $f_c$ >> B


texts: How to prevent Octave Jumping in ACF of Speech Signal?

texts:
I'm working with speech signals and my aim is to estimate the fundamental frequency
$\ F_0$   of this signal often called as "pitch".

The main idea is taking small blocks of the speech signal such that stationary can be assumed. Then calculating the autocorrelation function (ACF) of this block of a speech signal and finding index the global maximum of the ACF (except at zero) which refers to fundamental frequency.

But in the text it is stated that :

  The global maximum might not be at the lag corresponding to the true
  fundamental frequency but can possibly be an integer multiple of it.
  Due to this, the maximum can jump in consecutive frame between lags
  corresponding to multiples of T0 leading also to jumps in the
  F0-estimate. These effects are called octave-jumps.

My questions arise at this point: How does octave jumps occur? What is the possible reason? I know that ACF is a periodic function since original time sequence is periodic and in my opinion, this period equals to the block length of the original speech signal we are working on. When I investigate the interval, first period of the ACF, how can I decide whether the maximum is refers to the pitch or it is a maximum shifted from the consecutive period (block)? How can ı prevent from this effect?


texts: Identifying impulse based on experimental data

texts:
I have conducted an experiment in which the impulse was exerted on sound bowl (it was hit by a hammer). As a result I have obtained acceleration from accelerometer on 3 axis, calculated net acceleration ($\sqrt{x^2+y^2+z^2}$) and tried to obtain frequency components of the response.

What I need to do is to identify input impulse but I have no idea on how to do that.


texts: Problem in verifying zero blocking property of a zero of a multiple input multiple output system

texts:
I have a multiple input multiple output control system, and in particular a two inputs two outputs system. I want to find the zero state direction and the zero input direction of the system, which has a zero at $s=0$. So, I have foud them with the following code:

s = tf('s');

P = 1/(s+5);
C = 6/s;

S = 1/(1+P*C);
T = P*C/(1+P*C);

G_2 = [T S; S -S]; %transfer matrix

S = ss(G_2)
Smr = minreal(S)
A = Smr.A
B = Smr.B
C = Smr.C
D = Smr.D
trz = tzero(A,B,C,D)
RSM_1 = [0*eye(2)-A B; -C D]
rRSM_1 = rank(RSM_1)

null = null(RSM_1)

%verifying zero blocking property
x0 = [-0.4388; -0.1482];
u0 = [-0.0000; 0.8863];
t = linspace(0,5);
u = exp(t);
[y,x] = lsim(A,B*u0,C,D*u0,u,t,x0);
plot(t,y)
axis([0 5 -1 1])
xlabel('time, seconds')
title('response to x_{0} and u(t)=u_{0}e^{t}')

and I want to verify the zero blocking property, which means that if I apply an input at $u(t)=u_0e^{t}$ starting from $x(0)=x_0$, where $u_0$ is the input direction and $x_0$ is the state direction, the output is zero.

I have found it analytically in this way:

I know that the since the system zeros of a MIMO transfer function are those values $s = z$ which make the system matrix $S(z)$ loose rank, there exists at least a direction $u_0$ (zero input direction) and a n-dimensional vector $x_0$ (zero state direction), not simultaneously zero, such that

$\begin{pmatrix}
zI-A & -B\\
 C&D
\end{pmatrix}$ $\begin{bmatrix}
x_0\\
u_0
\end{bmatrix}$ $=$ $\begin{bmatrix}
0\\
0
\end{bmatrix}$

And by solving this system we obtain the zero state direction $x_0$ and the zero input direction $u_0$.

In order to make the question even clearer I can make an example using a different transfer function which has a simpler minimim realization. The example is taken from this : Here, and it can be found at the bottom of the paper.

The system in this case is:

$P(s)=\begin{bmatrix}
\frac{2}{s^2+3s+2} & \frac{2s}{s^2+3s+2}\\
 \frac{-2s}{s^2+3s+2}& \frac{-2}{s^2+3s+2}
\end{bmatrix}$

Its minimum realization is:

$A=\begin{bmatrix}
-1 &0  &0 \\
 0&-2 &0 \\
 0& 0 & -2
\end{bmatrix}$

$B=\begin{bmatrix}
2 &-2 \\
 -2& 4\\
 -4& 2
\end{bmatrix}$

$C=\begin{bmatrix}
1 &1  &0 \\
1 &0  &1
\end{bmatrix}$

$D=\begin{bmatrix}
0 &0 \\
 0& 0
\end{bmatrix}$

and the after using the same code, which can also be found in the paper, the simulation for verifying the zero blocking property is:

s = tf('s');
P = [2/(s^2+3*s+2) 2*s/(s^2+3*s+2); -2*s/(s^2+3*s+2) -2/(s^2+3*s+2)]
S = ss(P)
Smr = minreal(S)
A = Smr.A
B = Smr.B
C = Smr.C
D = Smr.D
trz = tzero(A,B,C,D)
RSM_1 = [eye(3)-A B; -C D]
rRSM_1 = rank(RSM_1)

null(RSM_1)

x0 = [0.5345; -0.5345; -0.5345];
u0 = [0.2673; -0.2673];
t = linspace(0,5);
u = exp(t);
[y,x] = lsim(A,B*u0,C,D*u0,u,t,x0);
plot(t,y),grid
axis([0 5 -1 1])
xlabel('time, seconds')
title('response to x_{0} and u(t)=u_{0}e^{t}')

and the output of this is:

which is correct, since the output coes to zero.

But, what I obtain with my code is:

I am still working on this, and by using the Matlab command isstable(T) it tells me that the closed loop is not stable. But if I run the step response, step(T) ,I find:

what am I doing wrong?


texts:
I normalized the spectrum of a time series and windowed the spectrum, then something strange happened. The code below can run in MATLAB or Octave.

figure;
n = 2e3;
t = 1:n;
m = n/2+1;
f = linspace(0,1,m);
x = rand(1,n) - 0.5;
y = fft(x);
subplot(2,2,1); plot(t,x,'k');
axis tight; xlabel('t/s'); title('original signal x');
subplot(2,2,2); plot(f,abs(y(1:m)),'k');
axis tight; xlabel('f/hz'); title('original spectrum');
y1 = y ./ abs(y);
k = n/10;
w = sin(linspace(0,pi/2,k));
y1(1:k) = y1(1:k) .* w;
y1(m:-1:m-k+1) = y1(m:-1:m-k+1) .* w;
y1(n-m+3:n) = y1(m-1:-1:2);
x1 = real(ifft(y1));
y2 = fft(x1);
x2 = real(ifft(y2));
subplot(2,2,3); plot(t,x1,'k',t,x2-x1,'r');
axis tight; xlabel('t/s'); title('whitened signal x1(black) and x2-x1(red)');
subplot(2,2,4); plot(f,abs(y1(1:m)),'k',f,abs(y2(1:m)),'r');
axis tight; xlabel('f/hz'); title('whitened spectrum y1(black) and y2(red)');

As shown in subplot 3, the red line is the difference between x1 and x2, meaning that they are exactly the same, but as shown in subplot 4, their spectrum (black and red, respectively) is different.
So, It's weird FFT of IFFT of a spectrum is different from the spectrum itself.


texts:
I read a paper about a brain-computer interface. In this paper the authors reported "each signal has been filtered with an 8-order band-pass Chebishev Type I filter which cut-off frequencies are 0.1 and 10 Hz and has been decimated according to the high cut-off frequency". I tried to design this filter with scipy:

import scipy.signal as signal
signal.cheby1(8,0.05,[0.1,10.0],btype='band',analog=0,output='ba')

The result was:

Warning: invalid value encountered in sqrt
(array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
    nan,  nan,  nan,  nan,  nan,  nan]), array([ nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,  nan,
    nan,  nan,  nan,  nan,  nan,  nan]))

I have no background in signal processing, so I actually don't know what I am doing. I don't know whether they used a IIR or FIR filter or whether I have to scale the cut-off frequencies or whether I'm using the wrong ripple. I hope you can help me.


texts: I'm currently working on a project to perform energy detection of an IEEE 802.15.4 using the block diagram illustrated below. However, I do have a little uncertainty about the algorithm. I'm wondering after finding the magnitude squared values at each FFT bin, do you then add each bin value up and divide by the number of bins (N) to get the average power? And also, do you just repeat this averaging M times? And why exactly do you repeat this averaging M times??


texts: Digital RF signal vector resize algorithm

texts:
I'm searching for an algorithm for resizing a RF sampled signal in order to draw it without losing any important maximum or minimum values.

The values of the vector are negative/positive around zero.


texts:
I want to print some line drawings on A3, but they aren't very large - around 600x800px.

I tried obvious functions in Gimp, like resize then enhancing with wavelet sharpener or unsharp mask, with selective gausian blur before it, but print results are no good.

They are not easy for vectorization either, but perhaps tweaking some tracing tool can give me better results then resizing bitmap. Before I start exploring this option, I thought to ask for advice:

Can anyone share a recipe, either by using programming algorithms (numpy/scipy or matlab) or by using tools available in Gimp or Photoshop, to accomplish quality resize of line drawing?


texts: How to distinguish the same object entering the camera field of view

texts:
We want to set photo cameras in the field and capture images of insects that visit flowers. We want to place these cameras above the flowers and then they should get triggered whenever an insect lands on the flower.
Is also important to find a way to distinguish if is the same insect just flying around and coming back in the next second or so (or a certain time frame). This is important because we would like to have an accurate estimation of how many unique insects visit a particular flower in a certain period of time. Sometimes an insect can create multiple false positives. That is, we can count that insect as multiple individuals, when is actually just the same individual.
I am wondering also how hard is to do this with conventional photo cameras that have incorporated motion sensors.
If this is not the right Stack Exchange site for such a question/discussion, please point me to a more suitable one.


texts: Why are some JPEG files of the same scene much larger than others?

texts:
I'm using a Foscam FI8910W ip camera to view a static scene under constant lighting conditions.  When I pull back a frame grab, it's about 35 KB in size.  I can do this over and over and it's always around 35 KB, but fluctuates somewhat due to various noises inherent to electronic image capture.  This random fluctuation is only of the order of 1 KB at most.

About every 2500 frames, the image size of a frame is suddenly in the order of 70 KB.  There's no gradual creep upwards if you're thinking thermal noise as the camera warms up.  1 frame will be 70 KB(ish) and then they return to the 35 KB sized frames.

This has happened before with another run looking at a different scene.  The common file size was 39 KB then, and out of 10,000 frames, 4 were of the order of 77 KB.  The image size histogram looked like this:-

Before you ask, I've managed to save one of these frames and it looks exactly like all the others bar the expected noise fluctuation.  They have roughly the same number of unique colours at around 23,000.  So it's not a moth randomly landing on the lens for exactly 1 frame and then flying off.  For completeness, I've done another run of images and this is an example typical image (the reflection is the IR illuminator):-

This is the anomaly image:-

You can see there is no difference.  Excuse the hippo.  I'm fairly familiar with the JPEG algorithm, and I can't see how this might happen other than a coding error with Foscam's implementation.  But, could there be something inherently chaotic within some of the JPEG transformation functions (like the  discrete cosine transform or  quantization)?  Statistically, one would expect a normal distribution of file size and this is what I see around 39 KB.  Then there are several outliers at 77 KB.  So it doesn't appear stochastic.

The reason this is in CS and not in hardware is I'm asking could this be a programming code phenomenon relating to the JPEG encoding algorithm? Seems unlikely, but the anomalies are random and infrequent and there is no human interaction with the device.  Is JPEG encoding stable?

The reason you might be unfamiliar with this phenomenon is that as the images look the same, no one really looks at the file sizes.  The file size is of critical importance to me so I noticed.   How can this happen roughly every 2500 frames?

Supplemental:-

Posting these images is not going to work guys, as the imgur software re-samples the uploaded files.  So whilst I posted 37K and 73K files, imgur has re-sampled both to 35K.  This seems to be a Stack Exchange problem which seems ironic for a  site dealing with image processing, data compression and analysis.

This is my processing of the images.  It's the normalised difference between a normal image and the anomaly.  The image is as you would expect, with JPEG noise in the high frequency regions.    This is an RGB image even though it looks monochrome.  There are 8000 unique colours in the colour cube (representing the noise).

Supplemental 2:-

As requested, 4 normal frames and 2 abnormal frames can the downloaded from sample frames.  It is a different scene, but the anomalous behaviour still occurred, so this proves that it's consistent.


texts: Can I study continuous time Fourier Transform and treat the rest as special cases

texts:
Say I learned the theoretical result of continuous time Fourier transform. And I want to extends some results(say "convolution rule") to

Lapace transform,
Z transform,
DTFT,
DFT,
Fourier sequence etc.

Then really, I just need to remember how one works and derive the rest. How to do this as a general principle? Is there any textbook that teaches this way? Thanks.


texts: How harmonics happens when you mix two signals?

texts:
Doing a search in quora about Why does Amplitude Modulation have sideband frequencies? i found an answer and within it, a user answered with the following statement:
"When you try mixing two signals, there will be something know as harmonics that will occur. This means that frqncy X (+) Frqncy Y will not just form frequency Z but some more frequencies that will be embedded into the system."
Can someone clarify and or exemplify the statement?


texts:
I wanted to know if there is a sound render like the renderers of 3d modeling software (vray ,  cycles etc).

I would like to experiment with sound in a controlled environment

if it does not exist which is what makes it difficult to create an audio renderer?


texts:
In Analog-to-Digital Converters(ADC), the signal is first sampled at a rate higher than or equal to the Nyquist, then quantized and encoded. In Analog-to-Informations (AIC), the sampling and compression are combined in one step, mathematically:
$$ y = \Phi \cdot x $$
where $y$ is compressed signal, $\Phi$ is the sensing matrix, and $x$ is the input signal. My questions are:
* Simulating the above equation using MATLAB is downsampling, First the input signal is read to obtain a vector of say $N$ samples representing the signal, and then multiplied by the matrix $\Phi$.

Does this mean that we sample the signal to obtain the vector $N$? If so, where is the benefit? How sampling and compression are combined in one step? How can the sampling and compression be combined to perform compressive sensing in a practical way?


texts: Generating Music Tones, How to handle when tones overlap? (16bit samples)

texts:
If i add the two tones then the result if two peaks are added together will be greater then the the sample max (2^15).

For the life of me I couldn't think of the right words to google. Think i have google block.


texts: Comparing two Accelerometer signals in time domain

texts:
As a research to answer my question, I've been reading several posts on this stack with a common thread of "comparing signals". The methods recurrently suggested were to use DTW, correlation or MSE etc.
The two signals I'm trying to compare are generated by piezoelectric single-axis accelerometer and they are periodic.
Out of 4 variables playing a role in signal amplitudes, two of the variables are different(A,B resp) and the rest are same(C,D resp).
My aim is to measure quantitatively how different the two signals are on the basis of variable B. At the moment, a qualitative way to go about this, according to my supervisors, is to look which among the two signals have the higher peaks.
I want to get a quantitative measure to compare the two signal's peaks or maybe the whole dataset which consists of 20,000 datapoints for a period of 10 sec. Are peaks/amplitudes the right way to go about this?


texts: Why are MATLAB and matplotlib spectrogram value different?

texts:
I'm doing a signal processing to put data into my neural network.

I noticed a difference between  MATLAB and  Python spectrogram Value. I don't know why and ask questions.

This is MATLAB code

[y, Fs] = audioread('left (1).wav'); %1sec 16000 sample read

nwin  = 256; %16ms
novlp = 128; %8ms
nfft  = 256;

left = spectrogram(y, nwin, novlp, nfft, Fs, 'yaxis');

256-fft_point, 256 hamming window, 50% over lab, 16000 sampling rate

And This is Python code

import numpy as np
import soundfile as sf
from matplotlib import mlab

audio_sample, sampling_rate = sf.read("left (1).wav")

win = np.hamming(256)
Sxx, f, t = mlab.specgram(x=audio_sample, NFFT=256, Fs=sampling_rate, window=win, noverlap=128, mode='complex',
                          sides='default', detrend='none')

np.savetxt('C:/Users/HOME/Desktop/left_spec (1).txt', Sxx, delimiter=',')

It is also 256-fft_point, 256 hamming window, 50% over lab, 16000 sampling rate

However, the values are totally different as shown below.

In addition to the first four values, all values show different values.

When I read the wav file, the value seems to be the problem of the spectrogram function because both MATLAB and Python are the same.

I also changed the parameters of Python spectrogram function, tried post-processing such as Sxx = 10 * np.log10 (Sxx), and tried the spectrogram function of the scipy kit, but not all the same values as MATLAB.

Please help me, thank you for reading.


texts:
So my question looks like this:

  Suppose that we have two LTI systems with impulse responses
  $$h_1(t) = \frac 12\delta(t-1)\quad\text{and}\quad h_2(t)=2\delta(t+1).$$
  Determine whether these systems are inverses of one another.

Here's what I've done but I'm stuck. I'd appreciate some help with this

$$h_1(t)\star h_2(t)\overset{?}{=}\delta(t)$$
\begin{align}
\int_{-\infty}^\infty h_1(\tau)h_2(t-\tau)d\tau &=\int_{-\infty}^\infty \frac12\delta(\tau-1)2\delta(t-\tau+1)d\tau\\
&=\int_{-\infty}^\infty \delta(\tau-1)\delta(t-\tau+1)d\tau\\
\end{align}

$$u = \tau - 1\implies \frac{du}{d\tau}=1\implies du = d\tau$$

$$\implies \int_{-\infty}^\infty h_1(\tau)h_2(t-\tau)d\tau = \int_{-\infty}^\infty \delta(u)\delta(t-u)d\tau$$


texts: Is there a way to separate 2 signals recorded on the same channel?

texts:
I have 2 low frequency signals on the same line that look like this:

I would like to separate them so that I obtain the two signals that can be seen here:

The sampling is not perfectly intercalated, so it's not as simple as taking one point each in alternation.

I looked at ICA but my understanding is that I would need 2 recordings of the 2 mixed signals in order to separate them.

Any ideas?


texts: Constant amplitude, uniform phase - what's the distribution of the complex signal then?

texts:
The well-known relationships for zero-mean circularly-symmetric complex Gaussian $z = a + jb = |z| \exp(j\varphi)$ signals are

the amplitudes $|z| = \sqrt{a^2 + b^2}$ are Rayleigh-distributed
the phases $\varphi = \arg(z)$ are uniformly distributed
the pdf of such signals is fully described by their complex covariance matrix $\mathbf{\Sigma}$, which is Wishart-distributed in case it is estimated from a number of independent samples ($\rightarrow \mathbf{\hat{\Sigma}}$)

My question is two-fold:

1) What's the distribution of $z$ in case the amplitudes are constant, e.g. all normalized to 1?

2) What's the distribution of $\mathbf{\hat{\Sigma}}$ in case the same normalization is applied, i.e. the amplitudes of its elements are removed?

I figure that it's then $z = \exp(j\varphi)$, with $\varphi \sim \mathcal{U}$, i.e. I am possibly looking for the distribution the exponential of a uniformly distributed variable follows. For $\mathbf{\Sigma}$, however, I have no idea.


texts: What is the equivalent smoothing function to running the same Gaussian 8 times?

texts: Suppose there is an image which is to be smoothed by convolving it with a Gaussian kernel with standard deviation $\sigma$. If the image is then smoothed with this kernel 8 consecutive times, is the total smoothing still Gaussian with some calculable standard deviation $\hat{\sigma}$? Or is it some other distribution instead?


texts: How does causality (i.e. unit step) affect the DTFT of a sine or cosine wave?

texts:
Tables of common Discrete-Time Fourier Transform pairs list the transform of a sine wave:

$ \sin(\omega_0\ n) $ and its transform:
$ -j\pi\ [d( \omega\ - \omega_0\ ) - d( \omega\ + \omega_0\ )] $

And the cosine:

$ \cos(\omega_0\ n) $ and its transform: $ \pi\ [d( \omega\ - \omega_0\ ) + d( \omega\ + \omega_0\ )] $

How might the results differ if the sin or cosine is causal? That is, I would like to determine the Fourier Transform of the following (sin or cosine multiplied by the unit step):

$ \sin(\omega_0\ n) u[n] $ or  $ \cos(\omega_0\ n) u[n] $

So regarding the DTFT of the signals above, is the only method to realize that multiplication in the time domain is convolution in the frequency domain? Or is there a more simple rule-of-thumb, or property, of the Fourier Transform that I am missing?

(This is not homework, but it is for studying, so it would be helpful to include resources along with answers.)


texts: Proper HW circuit to decode PWM signal using Arduino Mini

texts:
I want to properly connect a Radio Control Receiver signal to a Arduino Pro Mini, Pin 11 so I can detect the 1 ms to 2 ms PWM signal.  I calculate a value of 220 ohms will limit the circuit current to 22mA.  The Arduino strongly suggests using 20mA with a 40mA max on input pins.  Do I need a .1uF cap to ground?  Electrolytic cap?
 I would think adding a cap will cause the signal line to ride at a certain voltage level instead of having a sharp leading/trailing signal levels.
TIA


texts:
I am hoping that somebody might have a good idea for detecting low-frequency pulses in some rather noisy data (I'm working on improving the snr). I have repurposed a "Raspberry Shake 1D" to record the vibrations made by an adult filaria just below the skin and using a lowpass filter I was able to collect the following as a proof of concept to see if this was even possible. It seems to work.

https://raspberryshake.org/products/raspberry-shake-1d/

The issue here is I don't get any feedback from the device when it is attached and I need to download the stored data later as a "miniseed" file. The device is, however, able to send UDP packets in real-time and I would like to attach a display for a live picture of the noise filtered data. If possible I would like to add a pulse detection algorithm of some kind to make it more obvious when a distinct pulse of the appropriate frequency range, attack, and decay of this pulse is within some basic heuristic bounds.

I could take a deep learning approach but for that, I would need lots of tagged data that has not been collected yet. I am hoping for some form of fft like algorithm just to give me enough real-time feedback while placing the device to better collect my data.

Any thoughts?

Edit1:

What is being detected (shown in the graph) is a small subsonic vibration where the data leading into that waveform is the background ambient noise level, then the oscillations ramp up and back down, then the background level noise again as the graph ends. I think the frequency here is about 4Hz but that will likely vary, as well as the amplitude, so from the detection level, I may need something like what a PLL lock indicator circuit would do if this were a hardware problem, but then this is not a nice sine wave either. I have not done anything like this in software so I am not familiar with what software techniques might even be possible.

What I am envisioning is a python program receiving UDP packets, passing that data through a low-pass filter, then through some kind of FFT processing to find something near the center (4Hz?) frequency range, and the output from that to produce both my real-time graph and some indicator (noise or color change on the display) that the waveform's desired frequency and amplitude are present. It would need to ignore a signal over some specific amplitude, but I'll cross that bridge when I get something that mostly works. I have done this with the static recorded data, and visually identified the event, but I have no clue how to do this from live data in realtime.

btw - The note on "filaria" was only for context, so you could see the problem you are helping to solve, but it is not really important for this discussion. But since someone asked, it's an investigation of parasitic worm one acquires through being bitten by an infected insect, for which there is no diagnostic test in humans for the disease, no disease-specific symptoms, and no known cure for the disease. I'm trying to research a total enigma and these very minute vibrations are just one general property of the disease being investigated.


texts:
I would like to sum a filterbank made of bandpass butterworth filters to reconstruct an audio signal. Given group delays for each channel shown in the following plot . . .

How would I go about designing complementaty allpass filters to flatten this out somewhat? It doesn't have to be completely flat, just a bit flatter than it is. Sampling the impulse responses to make linear phase filters is not an option at the moment.

The Matlab code used to generate this is shown below. Note: I have multiplied the group delays by a factor of 2 as the signal is passed through the bank twice.

sr = 22050;

%Make octave wide filters
cfs = [125 250 500 1000 2000 4000];
bw = 1; %octaves

nBins = 2^12;
gd = zeros(nBins,numel(cfs));

for nn = 1:numel(cfs)
    cf = cfs(nn);
    loEdge = cf*(2^-bw/2);
    hiEdge = cf*(2^+bw/2);
    [b, a] = butter(2, [loEdge hiEdge]/(sr/2), 'bandpass');
    [gd(:,nn),f] = grpdelay(b,a,nBins,sr);
end

figure; semilogx(f,2*gd); xlim([10 8000])
legend(num2str(cfs')); xlabel ('F [Hz]'); ylabel('Group delay [samples]')


texts: How to evaluate expressions without explicitly computing a DFT?

texts:
I am trying to learn this part but I dont seem to understand.
if I have a finite length sequence, and a n-point DFT of it (in an interval)... is it possible to evaluate expressions of the  n-point DFT without computing them explicitly?

the exercise is:

Let $x[n]$ be finite-length sequence such that: $x[n] = \{2,1,1,0,3,2,0,3,4,6\}$ for $0<=n<=9$
    and $X[k]$ be the 10-point DFT of $x[n]$ for $0\le k\le9$. evaluate the following without explicitly computing the DFT:

$$I) \space  X[0]$$
$$II) \space \sum_{k=0}^9 X[k]$$
$$III) \space \sum_{k=0}^9  |X[k]|^2 $$
$$IV) \space \sum_{k=0}^9 X[k]e^{-j \frac {(4\pi)}{5}k}$$

and then it continues:

let $X[k]$ denote N-point DFT of an N-point sequence $x[n]$ ($n$ is even). Two N/2-point sequence are defined as:

$$g[n] = a_1  x[2n] + a_2  x[2n+1]  \space where \space  0\le n\le(N/2)-1$$
$$h[n] = a_3  x[2n] + a_4  x[2n+1]  \space where \space  0\le n\le(N/2)-1$$

constant $a_i$ satisfies: $a_1  a_4 \neq a_2  a_3$

let $G[k]$ and $H[k]$ denote the N/2-point DFT's of $g[n]$ and $h[n]$ respectively.
Determine N-point DFT $X[k]$ in terms of $G[k]$ and $H[k]$ .

any clues on the solution?

any input is highly appreciated.
thank you


texts: Matlab version rpmordermap and rpmfreqmap

texts: Are those functions only available in Matlab 2015b version or they are included in 2015a version also? Thanks in advance.


texts: Find Magnitude of a DFT signal using a Blackman window

texts:
Hi I've been given a signal made of a series of cosines. I have taken the DFT of the signal using a rectangular window (blue), hamming window (red) and a Blackman window (black). I have identified two of the cosine components of the signal. I was wondering how can I use the dB magnitude to find the actual magnitude of the frequency. This is so that i can recreate the signal by summing all of the cosines found.
I was also wondering why these frequencies can be seen on the blackman window but not the other 2 windows?
any help would be greatly appreciated!


texts: How can I determine the complex envelope of a signal when knowing the analytic signal?

texts:
Consider a signal $x(t) = \sum A_i \sin(2\pi(f_c+f_i) t + \theta_i),\, i=1\ldots N,$
and its analytic signal $Z_x(t) = x(t) + j\mathcal H[x(t)]$.

I want to describe the complex envelope (equivalent lowpass signal) as $x_{eq} (t) = r(t)\exp(j\Phi(t))$.

How do I know the amplitude $r(t)$ and the phase $\Phi(t)$ of the complex envelope?

As an example, let's consider: $x(t) = 1 sin(2\pi*11.5e9*t) + 2 sin(2\pi*11.6e9*t)$

In Matlab I would start the code with:

T = 1e-9;
f_s = 1e12;
t = [0 : 1/f_s : T]';
f_c = 11.5e9;

x_passband =  1*sin(2*pi*(f_c+0)*t) + 2*sin(2*pi*(f_c+0.1e9)*t);

x_analytic = hilbert(x_passband);

x_envelope = x_analytic .* exp(-1i*2*pif_ct);

The complex envelope should also be equal to:

amplitude = abs(x_analytic);

phase = angle(x_analytic);

x_envelope = amplitude * exp(j*phase);


texts: Discrepancy when calculating LTI system output using inverse z-Transform

texts:
I'm given a difference equation, $y[n]-0.4y[n-1]=x[n]$, and asked to find the natural response $y_n[n]$, forced response $y_f[n]$ and complete response $y[n]$ if $x[n]=4 (0.25)^nu[n]$ and $y[0]=0$.

By one approach I've read,
$$
\text{characteristic equation:}\quad z-0.4=0 \quad\therefore y_n[n]=A(0.4)^n \\
\\
y_f[n]=B(0.25)^n \quad\text{substitute back into diffeq:}\\
B(0.25)^n -0.4B(0.25)^{n-1}=B(0.25)^n - \frac{0.4}{0.25}B(0.25)^n= 4(0.25)^nu[n]\\
\therefore B=-\frac{20}{3}, \quad y_f[n]=-\frac{20}{3}(0.25)^nu[n]\\
\therefore y[n] = A(0.4)^n-\frac{20}{3}(0.25)^nu[n]\\
y[0]=0 \rightarrow A=\frac{20}{3} \\
y[n]= \frac{20}{3}(0.4)^n-\frac{20}{3}(0.25)^nu[n]\\
$$
Great! Alternatively, I should be able to determine $y[n]$ by taking the inverse z Transform of $Y(z)=H(z)X(z)$, which in this case is (I'm pretty sure)
$$
Y(z)= \frac{z}{z-0.4}\cdot 4 \frac{z}{z-0.25}\\
\begin{align}
\frac{Y(z)}{z} &= \frac{4z}{(z-0.4)(z-0.25)}=\frac{C}{z-0.4}+\frac{D}{z-0.25}\\
&=\frac{32}{3} \frac{1}{z-0.4} - \frac{20}{3}\frac{1}{z-0.25}\\
\end{align}\\
\therefore Y(z)=\frac{1}{3} \left[\frac{32z}{z-0.4} - \frac{20z}{z-0.25}\right]\\
\therefore y[n]= \left[ \frac{32}{3}(0.4)^n - \frac{20}{3}(0.25)^n \right] u[n]
$$
which is almost the same, but not quite.

Where did I go wrong? Also, in the second approach, how do you take into account the initial condition of $y[0]=0$? And thirdly, in the first approach, is the first term of $y[n]$ multiplied by $u[n]$ or not? If so, how do you show this?

(for anyone wondering, this is not homework for a course, just self study...)


texts: In what case do we need to increase the length of the signal compared to zero padding when using FFT?

texts: When do we have to increase the total length of the signal considered in time domain compared to just zero padding it when using FFT, this article http://www.bitweenie.com/listings/fft-zero-padding/ has a clear example and explains the difference. I understand the the concept of frequency bin and  FFT resolution but what is the waveform resolution when do we need to increase the length of signal in time domain?


texts: Use Scale Space Representation to Filter Single Image

texts:
Currently I hope to use scale space representation to filter one image. Features in one image can be filtered using an Gaussian smooth filter with one optimal sigma. It means different features in one image can be expressed best in different scale under scale space representation.

For example, I have one image with one tree in it. In the scale space representation, three sigma values are used and they are represented as sigma0, sigma1 and sigma2. The ground is best expressed in the smoothed image with sigma0 because it contains textures mainly. The branches are best expressed in the smoother image with sigma1 and the trunk is with the smoother image with sigma2. If I hope to filter the image, I hope that the filtered pixels for the group is from the smoothed image with sigma0.
The filtered pixels for the branches are from the smoothed image with sigma1. The filtered pixels for the trunk are from the smoothed image with sigma2.

It requires that I need to determine in which smoothed image one pixel is expressed best. Is this idea plausible?
I am trying to use differece-of-Gaussian of two successive smoothed images to perform the above task. Is there any other way to combine the three smoothed image?

I use MATLAB to implement the idea. The values of the three sigmas is 1.0, 2.0 and 3.0. The corresponding size of Gaussian kernel is 3, 5 and 7. I use the function fspecial to generate the kernel. Are the parameter reasonable? Please share your experience with the scale space representation to help me. You can provide some links to useful papers.


texts:
I have been stuck on this question for a while now. It has to do with vestigial sideband.

I wasn't sure if I should be dividing $H(\omega)$ graph values by 2 because only the positive side of the filter is given....

The answer that I am getting is:

We know that $VSB(t) = A_c f(t)\cos(\omega_c t) \ast h(t)$

then in frequency form I get:

$$VSB(\omega) = A_c\pi  0.5 (\delta(\omega + (\omega_c + \omega_m)) + \delta(\omega + (\omega_c - \omega_m)) + \delta(\omega - (\omega_c + \omega_m)) + \delta(\omega - (\omega_c - \omega_m))) H(\omega)$$

Now using the graph I am multiplying the appropriate $H(\omega)$ values:

$$VSB(\omega) = A_c \pi 0.5 (0.25 \delta(\omega + (\omega_c + \omega_m)) + 0.75 \delta(\omega + (\omega_c - \omega_m)) + 0.25 \delta(\omega - (\omega_c + \omega_m)) + 0.75 \delta(\omega - (\omega_c - \omega_m)))$$

where carrier wave is given by $A_c \cos(\omega_c t)$

I don't have the solution to this, so I have no idea if my answer or reasoning is correct...
It would be sooooo helpful if someone could confirm/clarify


texts: How should we characterize a non-LTI system? For example we have: $y[n]=x[3n]+x[2n]+x[n]$ which is clearly not LTI. Also, the impulse response will be $h[n]=3\delta[n]$ and if we take the DTFT of this system I believe $H(e^{j\omega})=3$ which is not very intuitive. How can we characterize this system?


texts: Difference between 2$\pi f$ and $\omega$ in Fourier transform

texts: What is the difference when we use $e^{-j2\pi f}$ and $e^{-j\omega n}$ for Fourier transformation?


texts:
For a project, I want to measure the contrast of a global image.

I want to compare different images to get the image with the best contrast.
I know about the Weber contrast, but I need a normalized contrast value (in the range between 0-1 or other).

Is there a simple method to do this in openCV?


texts: Hamilton-Jacobi-Bellman equation vs Riccati equation

texts: Most of the literature on Reinforcement Learning discuss Hamilton Jacobi Bellman equations for optimality. In dynamics however Ricatti equations are used. I am curious if there are any parallels between the two equations for optimality? How do one differ from the other?


texts:
I have some iq data recordings of a 2-FSK waveform that I've created periodograms from. How do you estimate the SNR of this waveform in the frequency domain? I've heard two opinions:
One:

Integrate the signal band power over the excursion
Find a clear area of the spectrum and integrate the noise power over the same bandwidth
Divide the two quantities: SNR = signal power/noise power
This seems weird because FSK is just two tones and integrating the signal over the excursion is just integrating a lot of noise and very little signal

Two:

Find the peak signal power of one of the tones
Estimate the median of the noise floor power
Divide peak signal by median of noise floor: SNR = peak / median
This also seems weird because depending on the snapshot of data, the amplitude of the tones is different. SNR isn't changing but this method doesn't give a constant SNR.

If I could collect the same iq data with the waveform turned off (just receiver noise), would it make the calculation any simpler?


texts: Filter design for frequencies in range 0-60Mhz

texts:
First of all excuse me for asking a seemingly vague question. But i am going nowhere in the problem i am trying to solve.

I am working on a powerline communication system where i need to upconvert in 28-56Mhz range and later downconvert it to 0-28Mhz so that it works in the coupler range.So i need to design a band pass filter in 28-56Mhz and a low pass filter in the 0-28Mhz range. I need to design an analog filter for the same. I have no experience of filter design. I have been looking online for resources but wherever i go there is an explanation on chebyshev, butterworth and  elliptic filters.Then you have 1st order,2nd order and higher order filters equations. My problem is from all this i am  unable to decide which type of filter to use and what order to use. So i need to understand practical design constraints and examples of how to design given certain constraints and then boiling down to a solution. So could anybody help me in how to approach the filter design problem and suggest either examples or resources to go forward.


texts: Side lobes not visible in FFT of rect function

texts:
I was making a plot to show side lobes for different spectral windows and I noticed if I take a simple fft of a rectangle window, I don't see any side lobes. I think I'm just making a simple mistake.

Here is the Python code I ran:

from matplotlib import pyplot as plt
import numpy as np
fAxis = np.linspace(-1/2,1/2,100)
rect = np.ones( 100 )
specWindow = abs(np.fft.fftshift(np.fft.fft(rect)))
plt.plot(specWindow)
plt.show()

Could anyone explain how to change my code slightly to see the side lobes I'd expect. I expect a discrete FT of a rectangle function to give me the Dirichlet kernel.


texts: Confusion about speech being piecewise stationary

texts:
From wikipedia:
"Modern general-purpose speech recognition systems are based on Hidden Markov Models. These are statistical models that output a sequence of symbols or quantities. HMMs are used in speech recognition because a speech signal can be viewed as a piecewise stationary signal or a short-time stationary signal. In a short time-scale (e.g., 10 milliseconds), speech can be approximated as a stationary process. Speech can be thought of as a Markov model for many stochastic purposes."

I get that if you break a speech signal into small enough parts, then a every part will make a very specific, basic sound, which I presume, has the same (or very similar) type of spectral envelope, as many other parts around it/close to it.. but I don't quite understand what the stationary process is, especially not in this context.

To be exact, I don't understand what exactly the stochastic process/it's random variables are, or even the motivation for describing a piece of signal as a stochastic process.
I mean the actual input is not a stochastic process, it's a concrete input, it's an observation - is the quote saying that basically any small enough piece of the input can be considered as a certain sample of some stationary process?

Do we call it stationary, because while we know that the small part has some sort of quality of being very repetitive, we don't know exactly where we've seperated it/where it starts? Similar to something like a stationary process of sin(t+Y), Y being a random variable.

I'm pretty confused about this whole thing- if I see a specific sample of white noise, is it correct to call it 'stationary'? I mean, the original distribution that generated this sample, that is a stationary process, but does that mean you call the sample/signal itself 'stationary'?

How exactly does the fact that sound is piecewise stationary imply it a good idea to use HMM? I simply get stuck at that - maybe a slow, thorough explanation of this implication (and what exactly it means that it's piecewise stationary) would be the most helpful to me.

Sorry if the questions are not understandable or not concrete enough. I'm simply very confused and don't know enough about this.


texts: How to send two packets in a point-to-point transmission with NOMA?

texts:
If a transmitter T1 has two packets M1 and M2 to send to a receiver R1, can T1 send M1 and M2 at the same time and using the same frequency based on Non-orthogonal multiple access (NOMA)?

If T1 sends x = P1*M1 + P2*M2 with allocated power P1 and P2. Then R1 will receive y1 = x*h1 + n1 where h1 is the channel gain between T1 and R1 and n1 is the noise. Can R1 decode both M1 and M2?


texts: Filtering and decimating input as sections in multiple stages

texts:
What would be my options to decimate a large input signal available only in sections (or blocks) due to DMA type of interface? The FIR should be about order of 18 only and the total decimation factor of 8. I've made a quick MATLAB function to illustrate this issue.

function decimFilt(varargin)
% Multi-Stage Block-By-Block Decimation Filtering

%% Check input args
enPlot = false;
if nargin>0
    enPlot = varargin{1};
end
%%

%% Definitions

% Sample rate and period
fs = 44.440e6;
ts = 1/fs;

% Input
N = 38400;              % Nbr of samples
B = 12;                 % Nbr of blocks
f1 = 0.0125*fs;
f2 = 0.25*fs;
n = (0:ts:N*ts-ts).';
x = sin(2*pi*f1*n) + ...
    sin(2*pi*f2*n+pi/6) + ...
    0.25*randn(size(n));
%%

%% Filter design
O = 18;
b = fir1(O,0.125);      % Linear-phase FIR filter.
%%

%% Filter and decimate entire signal in one block (REFERENCE).
m = 2;
%   1. stage
[yRef,Z1rf] = decFirFilt(b,x,m);
%   2. stage
[yRef,Z2rf] = decFirFilt(b,yRef,m);
%   2. stage
[yRef,Z3rf] = decFirFilt(b,yRef,m);
%%

%% Filtering and decimating the data as sections (BLOCK-BY-BLOCK)
xBlock = reshape(x,[N/B B]);
%   1. stage
y1m = zeros(N/(m*B),B);
Z1mf = zeros(O,B+1);        % Note! State vector needed per stage!
%   2. stage
y2m = zeros(N/(m*m*B),B);
Z2mf = zeros(O,B+1);        % Note! State vector needed per stage!
%   3. stage
y3m = zeros(N/(m*m*m*B),B);
Z3mf = zeros(O,B+1);        % Note! State vector needed per stage!
for i=1:B,
    % Load previous section final condition to the initial condition
    Zi = Z1mf(:,i);
    [y1m(:,i),Z1mf(:,i+1)] = decFirFiltCond(b,xBlock(:,i),m,Zi);
    % Load previous section final condition to the initial condition
    Zi = Z2mf(:,i);
    [y2m(:,i),Z2mf(:,i+1)] = decFirFiltCond(b,y1m(:,i),m,Zi);
    % Load previous section final condition to the initial condition
    Zi = Z3mf(:,i);
    [y3m(:,i),Z3mf(:,i+1)] = decFirFiltCond(b,y2m(:,i),m,Zi);
end
%%

%% Results
% Plot
if enPlot
figure, ...
    stem(yRef,'bo'), grid on, hold on, ...
    stem(y3m(:),'kx')
end

% MSE results
mse = mean(abs(y3m(:) - yRef).^2);
fprintf('decimFilt\n MSE:%.4e\n', mse)
%%

end

%% Internal functions
% Decimation filter
function [y, z] = decFirFilt(b, x, m)
[yUnDec, z] = filter(b,1,x);
y = yUnDec(1:m:end);
end

% Decimation filter with  access to initial and final conditions of the delays.
function [y, zo] = decFirFiltCond(b, x, m, zi)
[yUnDec, zo] = filter(b,1,x,zi);
y = yUnDec(1:m:end);
end
%%

I'm guessing some polyphase FIR or FFT-based overlap-add or -save implementations could also be used? What I've calculated the FFT-based overlap-add or -save would not be more efficient for this problem?


texts:
I have many information i want to be clear very well .

The defining characteristic of FIR filters optimal in the Chebyshev sense is that they minimize the maximum frequency-response error-magnitude over the frequency axis. In other terms, an optimal Chebyshev FIR filter is optimal in the minimax sense: The filter coefficients are chosen to minimize the worst-case error (maximum weighted error-magnitude ripple) over all frequencies.In another word, minimizing the worst case squared error induces a minimum Chebyshev error problem in some formulations.

In order to achieve an optimal equiripple complex FIR filter in the complex Chebyshev sense ,we can use an evolutionary algorithm.

What is the relation which allowed using the complex Chebyshev sense for optimisation equiripple complex FIR filter? What is effective evolutionary algorithm in this case?


texts: Red noise estimation in power spectra. How to model the red noise component of data and derive corresponding confidence levels.?

texts:
I have generated a power spectrum using fft (python) and now I want to estimate the confidence levels with respect to red noise. Something like this:
https://static-01.hindawi.com/articles/jcli/volume-2014/434960/figures/434960.fig.002a.jpg

The figure caption reads "Calculated spectrum (solid black), red noise curve (green solid), and the curves indicating the upper (red dashed) and lower (blue dashed) bounds of 95% confidence level, for WCA (a) and ECA (b) subregions. The red noise spectrum is computed using the Markov processes".

I searched the internet for the last two days with keywords such as "statistical significance of fourier spectrum" and "red noise in fourier spectra", but without much success. Can some please tell me the steps to derive the confidence levels for a power spectrum?


texts: can MFCC coefficients alone be used to do voice detection?

texts:
Aim of my project is to differentiate between cases of silence, single active speaker and multiple ( simultaneous ) active speakers using as basic ( but effective ) logic as possible. Is it possible to achieve this without use of extensive machine learning algorithms ??

Actually, i don't need to do speaker "identification". An audio stream will be taken as an input and I have to just say whether that audio consists of a single speaker, multiple speakers or silence... I have to implement this in C (real time), hence machine learning (unsupervised) would increase the computational complexity.


texts:
I'm trying to find the equation of a shelving filter of the second order. I easily found the equation for a first order filter on wikipedia with a transition region which have a 6db per octave slope. The equations are :

coef^2 * (1 + (f / (coef * f0))^2) / (1 + (f/f0)^2) for the lowshelf
(1 + (f / (coef * f0))^2) / (1 + (f/f0)^2) for the highshelf

With :

coef = multiplicator coef (gain = 10 log10 (coef))

f = frequency

f0 = cutoff frequency

But now, how can I use a kind of Q factor to modify the slope of the transition region ?

Of course, as I'm using already a frequency signal and not a temporal signal, I can't use the equations made for the z transform...

Thank you for all your answers


texts:
I'm a beginner in DSP and I'm going through the textbook of Oppenheim's Discrete Time Signal Processing. There are two figures in the text, one which I can visualize, and the other I can't. The first figure shows the polyphase decomposition of an impulse response function:

I get why the $z^{-n}$ is placed after the upsampling operation. After the zeros are added during the upsampling operation, appropriate delays need to be added so that the sum gives the original impulse response. The next step taken in the textbook isn't clear to me. I do not get why the $z$ delays are added before the polyphase components in the second figure, as shown below:

Doing so suggests the delay operation is commutative with the expansion, but I cannot convince myself that it is right. I tried doing that for the first figure and I do not get the same results. Is there something I'm misunderstanding here?

Thanks.


texts: Mathematical relationship between phase modulated (PM) signal and the complex phase/magnitude of the frequency component at the carrier frequency?

texts:
What is the mathematical relationship between phase modulated (PM) signal and the complex phase/magnitude of the frequency component at the carrier frequency used?
This is partially related to the question that I asked on phase noise, but it seems that the complex magnitude of this component is the same for all PM signals using that carrier, just that the phase is rotated. Is it something like the average instantaneous phase difference?


texts: Z Transform - Do i always need 2 poles for every "peak"?

texts:
I am quite new to digital signal processing (and also Z transforms). I am reading about frequency response modelling and have some questions

do we always need a pair of poles for each "peak" in the frequency response ?
why is this so ?
what happens if i dont have a pair of pole for each peak ?


texts: Confusion regarding frequency of Discrete Time Periodic Signal

texts:
While Analyzing the DFT Plot for a signal with N samples, say we find a peak in the magnitude of the DFT plot at index $k$, this implies that our signal has a high amount of similarity with an exponential phasor rotating at $\omega = (2\pi/N)k$ or equivalently a period of $N/k$ samples.
But since we are in discrete time, the fundamental period must be an integer i.e. let the fundamental period be $T$ samples
$ \implies T/m = N/k$ (where T and m are relatively prime)
Hence shouldn't the period be $Nm/k$ and hence the frequency would be $(2\pi/N)(k/m)$  ? (where $m$ is the smallest integer that makes $T$ a positive integer.)
Thank you


texts: Equivalent of the Shepp-Logan phantom for video denoising

texts:
The Shepp-Logan phantom is a standard test image used in reconstruction algorithms as well as image denoising.

Does an equivalent image sequence exist for testing video denoising algorithms? I know that there are plenty of "natural" videos out there (one of two table tennis players springs to mind), but I'm looking for a more artificial test sequence.


texts:
Plotting the FFT of a signal with Matlab e.g. by the code:

y_fft = abs(fft(y));
f = fs*(0:length(y)-1)/length(y);
plot(f,y_fft);

where y is the discrete signal and 'fs' is sample rate. creates a plot where the curve can be mirrored in the middle of the x-axis e.g.

Why can the FFT always be mirrored in the middle of the x axis?


texts: Pole Magnitude and Damping Ratio relationship

texts: I know that the damping ratio of a system is defined by the angle of the pole, calculated with respect to the left hand side $x$-axis. Could one infer though, that if the magnitude of the poles is small (i.e. the conjugate poles lie closer to the origin), then the damping ratio of the system would be small ,and result in a slower system response?


texts:
I'm seeing something very curious in the Matlab code below.

X = fft(x);
Xt = X;
Xt( (end/2+1):end ) = 0;
xr = 2 * real( ifft( Xt ) );
Xr = fft(xr);

where x is a real zero mean sinus function sampled at Fs Hertz according to Nyquist.

The spectra X and Xr should be the same. However, Xr(Fs/2) is not equal to X(Fs/2). Also, xr and x are slightly deviating.

Am I missing somehting here?


texts:
I have to find Fourier coefficients of $x(t) = |A \cos (2 \pi f t )|$. The problem also give me $y(t) = |x(t)|$. To find Fourier coefficients I wrote

$$ \frac{1}{T_0} \int_{0}^{T_0} |A \cos (2 \pi f t )| e^{2 \pi k f t}$$

Solving this I found that $k$ should be equal to  $1$ or $-1$. First I assume that $=1$ and I found the result $A/2 $ but on the book should be $2A/3\pi$. $$ T=\frac{T_0}{2} $$ Can someone help me?


texts: Result of conjugate symmetry property of DFT

texts:
I know one of the properties of DFT for real-valued time series is conjugate symmetry.

But what does it imply? In the textbook it says that for a DFT of the length M, this makes M/2-1 spectral coefficients redundant. Could someone please explain this to me? How do we conclude this?


texts: understanding FIR SRRC filter with memory implementation in C

texts:
I translated the Matlab QPSK receiver/transmitter into real-time C code that I have running on a PlutoSDR platform. QPSK receiver/transmitter
I wrote the square root raised cosine filter for the TX chain and it worked without any issues. I was able to successfully encode my message data and continuously transmit it then receive it and successfully decode it. The code for the root raised cosine filter in the TX chain is below:
static const double filterCoeff[21] = { -0.00045472051487622253, 0.0035368955557498576,
-0.0071456080909122616, 0.0075790619051782785, 0.0021436824272736714,
-0.010610686667249576, 0.0300115539818315, -0.053053433336247971,
-0.075028884954578726, 0.40916871463405252, 0.80373860039798006,
0.40916871463405252, -0.075028884954578726, -0.053053433336247971,
0.0300115539818315, -0.010610686667249576, 0.0021436824272736714,
0.0075790619051782785, -0.0071456080909122616, 0.0035368955557498576,
-0.00045472051487622253 };

void use_rrc_filter(int taps, int siglen, double complex* input, double complex* output, double*
rrc_filter_taps)
{
   double complex sum = 0;

   for(int i = 0; i < siglen; i++)
   {
       for(int j = 0; j <= i && j < taps; j++)
       {
           sum += input[i-j]*rrc_filter_taps[j];
       }

       output[i] = sum;
       sum = 0 + 0*I;
   }
}

The number of taps is given by: 10*2 + 1 = 21, where 10 is the filter span in symbols, 2 is the number of samples per symbol. The filter taps were generated by Matlab.
The problem arose when I tried to use this same code for the filter in the RX chain. It did not work. I realized I needed a filter that has "memory" and takes the previously received values into account (non causal). I was able to get MATLAB to generate C code for this filter which worked for me. The code is below:
typedef struct rcRXfilter
{
   int PhaseIdx;
   int CoeffIdx;
   int TapDelayIndex;
   double complex Sum;
   double complex StatesBuffer[RRC_TAPS-1];
   double complex output[FRAME_SIZE*2];
   bool isInitialized;
}rcRXfilter;

void initialize_rrc_rx_filt(rcRXfilter* obj)
{
  if (obj->isInitialized != 1)
  {
      obj->CoeffIdx = 0;
      obj->PhaseIdx = 0;
      obj->TapDelayIndex = 0;
      obj->Sum = 0.0;
      for (int i = 0; i < 20; i++)
      {
        obj->StatesBuffer[i] = 0.0;
      }

      obj->isInitialized = 1;
  }
}

void rrc_filt(double complex* input, double complex* output, rcRXfilter* obj)
{
   if (obj->isInitialized != 1)
   {
       initialize_rrc_rx_filt(obj);
   }

   int i = 0;
   int iIdx = 0;
   int jIdx = 0;
   int outBufIdx = 0;
   int curTapIdx = obj->TapDelayIndex;
   int phaseIdx = obj->PhaseIdx;
   int cffIdx = obj->CoeffIdx;
   int maxWindow = (phaseIdx + 1) * 20;

   for (iIdx = 0; iIdx < FRAME_SIZE*2; iIdx++)
   {
      obj->Sum += input[i] * filterCoeff[cffIdx];
      cffIdx++;

   for (jIdx = curTapIdx + 1; jIdx < maxWindow; jIdx++)
   {
      obj->Sum += obj->StatesBuffer[jIdx] * filterCoeff[cffIdx];
      cffIdx++;
   }

   for (jIdx = maxWindow - 20; jIdx <= curTapIdx; jIdx++)
   {
      obj->Sum += obj->StatesBuffer[jIdx] * filterCoeff[cffIdx];
      cffIdx++;
   }

   obj->StatesBuffer[curTapIdx] = input[i];
   i++;
   curTapIdx += 20;

   if (curTapIdx >= 20)
   {
      curTapIdx -= 20;
   }

   phaseIdx++;

   if (phaseIdx < 1)
   {
      maxWindow += 20;
   }
   else
   {
      obj->output[outBufIdx] = obj->Sum;
      outBufIdx++;
      obj->Sum = 0.0;
      phaseIdx = 0;
      cffIdx = 0;
      curTapIdx--;

      if (curTapIdx < 0)
      {
        curTapIdx += 20;
      }
      maxWindow = 20;
   }
 }

  obj->TapDelayIndex = curTapIdx;
  obj->CoeffIdx = cffIdx;
  obj->PhaseIdx = phaseIdx;

  for (i = 0; i < FRAME_SIZE*2; i++)
  {
    output[i] = obj->output[i];
  }
}

Can someone help me understand this code? I know it is performing downsampling but I'm having trouble understanding exactly what it is doing. Any equations illustrating what is happening and/or a block diagram would be ideal. Thank you!


texts: Generating a Sampled Square Wave in Matlab

texts:
I'm a newbie on use of MATLAB and also in signals and systems. I need clarification and guide (with fundamental explanations please) on above.

I want to generate a sampled square wave in MATLAB with following criteria;

Signal duration $T=5\mathrm{s}$, Sample frequency $f_s=40\mathrm{Hz}$. there are supposed to be $f_sT$ signal points (what does this mean please).
$\tau=1/100$. the square wave signals are zeros for the first $\frac{f_sT(1-\tau)}{2}$ points and also for the last $\frac{f_sT(1-\tau)}{2}$ points. they are one otherwise.

What I have done so far:

clear all

fs = 40;
t = [0:1/fs:5];
dc = 50;
x = square(2*pi*t);
%f=1 for a square wave

plot (t,x);

Am I right? I have no idea how to use tau, and what signal points imply here. Kindly assist. Thank you.


texts: Discrepancy between Gaussian FFT and its Fourier Transform

texts: I am trying to do the FFT of a Gaussian signal and comparing it to the theoretical Fourier transform. For infinitely small time step $dt$ and infinitely long signal length $T$, the 2 should become identical right?


texts:
I have a motion signal that you can think of as a sequence of pairs (x,y); what would be the best approach to decimate this sequence at real time?
I haven't found anything similar on the web and I'm still learning about dsp.

Thanks


texts: Is there a convolution mistake in my method?

texts:
I have the input signal $x(t)$

And impulse response $h(t)=20 e^{-1000t} u(t)$ in which u(t) is the unit step function.

When I try a convolution, I thought the solutions would be something like:

$
\begin{array}{ll}
  \int\limits_{0}^{t} 2\cdot 20e^{-1000(t-\tau)}d\tau = \frac{1}{25}(1-e^{-1000t}),&  0 \le t <  2 \\
  \int\limits_{0}^{2} 2\cdot 20e^{-1000(t-\tau)}d\tau - \int\limits_{2}^{t}  20e^{-1000(t-\tau)}d\tau = \frac{3}{50}e^{2000-1000t}-\frac{e^{-1000t}}{25}-\frac{1}{50},& 2\le t <  3 \\
 \int\limits_{0}^{2} 2\cdot 20e^{-1000(t-\tau)}d\tau - \int\limits_{2}^{3}  20e^{-1000(t-\tau)}d\tau = \frac{3}{50}e^{2000-1000t}-\frac{1}{50}e^{3000-1000t}-\frac{e^{-1000t}}{25},&  t \geq 3 \\
 \end{array}
$

But this did  not look like the answer in the solution manual, so did I do something wrong here?

Answer in solution manual:

Alternative solution in solution manual:


texts:
I'm currently trying to decide what classifier to use for my snore detection algorithm. Which uses some features like Zero cross rating, frequency band energy, first formants...

The classifiers which have been pointed to are:

Support Vector Machine
Naive Bayes Classifier
Knn Nearest Neighbour
Linear or Quadratic Classifier

However, could somebody describe what situation each classifier is used in,  or point me to a good paper, website for article which describes each classifier well.. I find the Mathworks website hard to understand the difference between each.

thanks


texts: Impulse response (general form for linear systems as two-variables function $h(t,\tau)$) applied to Time-invariant systems

texts:
If a system is linear, then the operator $S$ mapping input signals into output signals - i.e. $y(t)=S\{u(t)\}(t)$ - is the integral of the input weighted by the impulse response:
$$ y(t) = \int_{-\infty}^{+\infty} h(t,\tau) u(\tau) \mathrm{d}\tau $$
where $h(t,\tau):=S\{\delta(t-\tau)\}$, response of the system to an impulse in $\tau$, is a function of two variables: the second one $\tau$ refers to the time at which the input is applied while the first one $t$ refers to the time at which the output is observed.

If this system is also time-invariant, the system operator commutes with the time shift operator:
$$S\{u(t-T)\}(t) = S\{u(t)\}(t-T) $$
that is, time shifting the input, a time shifted version of the output is produced.
Linear time invariant systems have an impulse response function which depends upon a single variable, the difference between $t$ and $\tau$:
$$h(t,\tau)=\tilde h(t-\tau)$$ where the $\sim$ symbol is used just to underline that $h$ and $\tilde h$ cannot be the same, because they have a different number of arguments.
How can I show that the impulse response can be expressed as a function depending only on $t-\tau$ starting from the definition of time invariance?


texts:
MUSIC algorithm has been known to provide super-resolution imaging capabilities, and it is implemented in target recognition radar system, however, I couldn't find a definite expression for its resolution (minimum spacing between two components).

is there a known expression for the resolution improvement provided by MUSIC algorithm over -for example- the FFT method?


texts: What's the essential bandwidth of the unit step function?

texts: The Fourier spectrum is in the Figure, how to find the essential bandwidth?


texts: Understanding and Correcting phase shift in General Linear Phase FIR filters

texts:
I have been learning about General Linear Phase Filters and the four standard FIR implementations.

I was wondering in the case $\beta = \frac{\pi}{2}$ or $\beta = \frac{3\pi}{2}$ if there was a way to correct the phase distortion caused by the filter? More generally, what does the the $\beta$ term actually do when looking at the signal in the time domain? It has the potential to multiply in a complex constant containing an imaginary term but I don't understand what that means.


texts:
FSK modulation encodes information into the frequency of the signals. But when I checked the reference books, the formal expression for a frequency modulated signal is:
$x_i(t) = \alpha cos(2\pi f_it + \theta_i)$
Which means, for a binary FSK, we can consider $x_1(t)$ and $x_2(t)$ with different frequencies $f_1$ and $f_2$ (which make sense) and different phase $\theta_1$ and $\theta_2$ (which is odd to me), respectively.
(1) My first question is, why should we allow the signals have different phases? I cant see a reason for this and as far as I know, different phases result in phase discontinuity of the transmit signal.
(2) And the second question is, how this difference in phase is related to coherent and non-coherent receiver?


texts:
I'm new here and I hope you can help me.
I started studying the Fourier Transform now at University and I have a lot of doubts about this subject.
My teacher sent to us one list of exercises and I don't know how I can solve the Fourier Transform of this function:

$$
x(t) = \delta \left(\frac{t-T}{a}\right)
$$

I don't know how to use the properties of scaling and shifting in this case, it would be great if you could explain to me the steps to solve this question.
Thanks for your help :)


texts: Standard deviation of Gaussian corresponding to perfect low pass filter [1,1]

texts:
For a time discrete signal, the rect-filter (for 1D signal [1,1]) is a perfect low pass filter. I was wondering, what is the best gaussian approximation of this filter?

For a current problem, I cannot use binomial filters but have to use Gaussian filters for smoothing. Ths question is thus, what is the standard deviation of a Gaussian filter best approximating the [1,1] rect-filter. 'Best' could e.g. mean is has similar probperties in the frequency domain.


texts: Windowing before viewing a filter frequency response?

texts:
It is well known that it is a good idea to apply windowing to a signal, before obtaining its DFT and viewing its frequency response.

But is it a good idea to apply windowing to the impulse response of a filter, before viewing its frequency response?
Does the answer to this question depend on whether the filter is FIR or IIR?

I could perhaps see why windowing could be better than just truncating an IIR filter. But I can't decide whether the same applies to an FIR filter.

However, even in the case of an IIR filter, we usually don't do this. Instead we usually just plot the response as a function of the IIR filter's transfer function numerator & denominator coefficients (e.g. using freqz in Python or Matlab).

So is it ever a good idea to window the impulse response of a filter?


texts: How to determine stopband of discrete Gaussian, stdev sigma, support N

texts:
I would like to analyse a gaussian smoothing kernel with a set standard deviation and support (let's say, in MATLAB, fspecial('gaussian', [5 1], 1.3) so sigma is 1.3 and support is 5) in the DTFT domain to determine the stopband.

My understanding is that the stopband is when the magnitude of the impulse response is below -3 dB. However, I don't have a clue how to determine this given a sigma and support value.

I'm particularly interested to see if the DTFT quantifies the artifact due to the compact support as I would like to understand this a little better.


texts: Power spectrum of 2D image: result interpretation

texts:
Let's say I have an image, which is NxN pixels in size. If I calculate the fourier power spectrum of this image, I get NxN values, with the highest distinguishable frequency at +- N/2 in each direction.

Let's also say I wanted to get a histogram of power over frequency, because I want to know which frequency carries which amount of power.
Then I'd create concentric circles, with the origin of the fft. Bigger radius (and bigger distance from the origin), higher frequency, etc.

But now the questions start:

If the input image is not a square, but a rectangle, the highest frequency still has to be N/2 and M/2, respectively, and again both + and -.

In my opinion, it does not make any sense to still use circles for calculating the histogram. I'd have to use ellipses, wouldn't I?

Imagine a circle/ellipse now, going in each direction to the max frequency.
Then there are four areas (area of the square minus the area of the circle) in each corner, from which I don't get any contribution to this histogram, because the absolute distance is larger than my max frequency.

Where's the error in that?

And last, I read several books and articles in the past few weeks. I read a passage about the same problem, but they didn't use circles/ellipses, they used a crooked square, i.e. a rhombus, with each corner on one of the axis.
The problem is, I can neither find the passage again nor remember the actual text (for full text search).
Does this approach make any sense? What about the four triangle-shaped areas in this case?

Regards and thanks for your help.


texts: How to Realize the Sigma Point Sampling Function in Unscented Kalman Filter?

texts:
Recently I'm learning the unscented kalman filter (UKF). When designing the unscented kalman filter, it involves a non-linear function to generate the sigma points and then use the system non-linear function to transform those sigma points.
However, I'm still confused about why we need to select sigma points to compute the estimate value and how the sigma points selecting function affects the efficiency of UKF.

Here is the pdf book I'm learning now:

https://drive.google.com/file/d/0By_SW19c1BfhSVFzNHc0SjduNzg/view

In this book, it mentions the sigma function in Chapter 10.3


texts: Are capacity and spectral efficiency identical for a dicrete-time digital signal?

texts:
I saw this answer and saw on the wikipedia page that both spectral efficiency and capacity can be measured in bpcu i.e. 'bits/symbol'. I just want to confirm then that they are identical for a discrete-time digital signal.

This diagram is plotting $\log_2(1+SNR)$ but referring to it as 'capacity', which confirms my suspicion because this should typically be spectral efficiency before being multiplied with bandwidth to get capacity. I consistently see diagrams showing 'spectral efficiencies' and 'capacities' showing this range of values.
I do not know where the $1/2$ factor came from in the answer I linked, or why they've used variance of the discrete noise in the SNR, but it appears that capacity and spectral efficiency are the same for a digital signal.
My guess would be that 'channel use' means per sample. Imagine an OFDM with 64 subcarriers using 16QAM and pretend none of the carriers are used for pilots or guard bands. The total number of bits per symbol will be 4*64 i.e. 256, but because there are 64 samples in a symbol and hence 64 frequencies used, if you divide that by 64 you get 4. Therefore spectral efficiency = channel use in this scenario?
If 48 of the 64 subcarriers are used to transmit data and the rest are pilots and guard bands, then I'm guessing that the spectral efficiency changes to 192/64 and changes to 192/80 if you include the cyclic prefix i.e. 2.4 bpcu before taking into account error coding of 3/4 where it becomes 1.8bpcu. And modulation efficiency is 4 bpcu now? i.e. 192/48.
Is this because there is no concept of a capacity as distinct from spectral efficiency in a discrete-time signal?
If you calculate the spectral efficiency of the OFDM signal using the bps/Hz method, you get the same result i.e. $\log_2\left(1+\frac{S}{N}\right)$, but the capacity is instead $B\log_2\left(1+\frac{S}{N}\right)$ i.e. 20MHz x 4 bits = 80Mbps. You can also calculate this capacity if you know the capacity in bpcu and know the time interval of the channel use, then you do 4*sample freuquency = 4 x 20MHz = 80Mbps


texts: Discrete-time, second-order section, direct-form II filter

texts:
I found dfilt.df2sos in Matlab and trying to figure it out how do I do that in Python? Is there any function that works similarly with the Matlab function?
I also tried to browse but I could not find a similar function in Scipy.

More dfilt.df2sos https://www.mathworks.com/help/signal/ref/dfilt.df2sos.html

Best,