texts: Error in using Kalman Filter for 2D Position Estimation in Python

texts:
This is my first question on DSP Stack exchange, so I apologise if it is poorly worded. I have some positioning data from a vehicle (GPX Format, collected through Strava) and want to use a Kalman filter as a first step to track the position of the car and determine the accuracy of the GPS measurement. However, the output estimate of the Kalman filter from my implementation seems completely wrong and absolutely does not match the position information and plot from the data. Can someone help me figure out what is wrong with this implementation and provide some methods to fix this problem?
Appreciate all the help in advance. Thanks!

As a reference, I used the code given in the below link and correspondingly modified it based on my requirements: https://stackoverflow.com/questions/13901997/kalman-2d-filter-in-pythonn.
Background: I only have a basic understanding of the working of the Kalman filter and am a new user to Python, but for this implementation, I've considered a constant velocity model with states as Position and Velocity, time step is assumed to be 1 (Considering GPS updates at 1Hz), the measurement matrix only considers the position information and the actual measurement gives the corresponding longitude and latitude values.
The test GPX file was obtained from the following link: https://github.com/stevenvandorpe/testdata/blob/master/gps_coordinates/gpx/my_run_001.gpx

My implementation in Python:

import gpxpy
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

with open('test3.gpx') as fh:
    gpx_file = gpxpy.parse(fh)
segment = gpx_file.tracks[0].segments[0]
coords = pd.DataFrame([
    {'lat': p.latitude,
     'lon': p.longitude,
     'ele': p.elevation,
     'time': p.time} for p in segment.points])
coords.head(3)
plt.plot(coords.lon[::36], coords.lat[::36], 'ro')
plt.show()

def kalman_xy(x, P, measurement, R,
              Q = np.array(np.eye(4))):

    return kalman(x, P, measurement, R, Q,
                  F=np.array([[1.0, 0.0, 1.0, 0.0],
                              [0.0, 1.0, 0.0, 1.0],
                              [0.0, 0.0, 1.0, 0.0],
                              [0.0, 0.0, 0.0, 1.0]]),
                  H=np.array([[1.0, 0.0, 0.0, 0.0],
                              [0.0, 1.0, 0.0, 0.0]]))

def kalman(x, P, measurement, R, Q, F, H):

    y = np.array(measurement).T - np.dot(H,x)
    S = H.dot(P).dot(H.T) + R  # residual convariance
    K = np.dot((P.dot(H.T)), np.linalg.pinv(S))
    x = x + K.dot(y)
    I = np.array(np.eye(F.shape[0]))  # identity matrix
    P = np.dot((I - np.dot(K,H)),P)

    # PREDICT x, P
    x = np.dot(F,x)
    P = F.dot(P).dot(F.T) + Q

    return x, P

def demo_kalman_xy():
    x = np.array([[100.0, 0.0, 0.0, 100.0]]).T
    P = np.array(np.eye(4))*1000 # initial uncertainty
    plt.plot(coords.lon[::36], coords.lat[::36], 'ro')
    result = []
    R = 0.01**2
    for meas in zip(coords.lon, coords.lat):
        x, P = kalman_xy(x, P, meas, R)
        result.append((x[:2]).tolist())
    kalman_x, kalman_y = zip(*result)
    plt.plot(kalman_x, kalman_y, 'g-')
    plt.show()

demo_kalman_xy()


texts: Determining asymptotic stability using transfer function?

texts:
In an exam task, I am asked to determine the transfer function of the following direct-time system and decide whether it's stable.

I think this system is canonical and the amplifiers 'on top' correspond to $b_2$, $b_1$ and $b_0$ from left to right in the following formula for the transfer function. Also, the amplifiers 'on the bottom' correspond to $-a_2$ and $-a_1$ similarly from left to right:

$$ H(z) = \frac{b_0 + b_1 z^{-1} + b_2 z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}} $$

Therefore the transfer function of this system is

$$ H(z) = \frac{z^{-1} + 0.4 z^{-2}}{1 + (-0.5) z^{-1} + 0.06 z^{-2}} = \frac{z + 0.4}{z^2 + (-0.5)z + 0.06} $$

So it has a single zero at $z = -0.4$ and its poles are $p_1 = 0.3$ and $p_2 = 0.2$

I know I can draw the pole-zero plot for this function that will look like this:

How do I deduce the stability of the system from here?

I have learned things before like given the eigenvalues $\lambda_i$ of the system's $\underline{\underline{A}}$ matrix, a discrete-time system is asymptotically stable if $\forall \lambda_i : |\lambda_i| < 1$. Similarly, in continuous time, the criterion is $\forall \lambda_i : \mathrm{Re}\{\lambda_i\} < 0$.

But now I have poles, not eigenvalues. I have only found that if there is a pole that is to the right of the imaginary axis on the pole-zero plot, the system is unstable. However, here, I think this system should be asymptotaically stable (therefore also BIBO stable), but it has a pole in such a place.

Why is the transfer function given this way for such systems?


texts:
This is a problem I encountered some time ago while doing my thesis, which still now, I am not sure to understand.
By then, I was implementing the circular correlation in an FPGA using this Xilinx FFT IP core as the acquisition block for a GNSS receiver. Aiming at reducing the required FPGA resources, I studied the impact of some of the IP configurations (scaling and truncation), and of different input signal amplitudes (but with a fixed data width of 8 bits) on the correlation in terms of distortion and SNR reduction.
I computed the FFT of a clean GPS C/A code using the Xilinx IP core, and then I computed its auto-correlation (ACF) with Matlab (i.e. $ACF=|IFFT_{Matlab}\left(FFT_{Xilinx} \cdot  FFT_{Xilinx}^* \right)|$). Similarly, I computed the FFT of a different GPS C/A code using the Xilinx IP core, and then I computed the cross-correlation (XCF) with the previous code, also using Matlab .
The impact on the auto-correlation (Fig. 8.12) was what I expected, an increase of the "noise floor" caused (I believe) by the accumulated quantization noise generated after each butterfly stage of the FFT. However, for two IP configurations, the cross-correlation also showed a peak (Figs. 8.12c and 8.12d). Although it's true that these cross-correlation peaks are very small  (probably well below the thermal noise of a real signal), a non present GPS satellite could be flag as present if the SNR is large enough, or if a lot incoherent averaging is performed.
So my question is, can the truncation/rounding in a fixed-point FFT cause a deterministic signal that is common in signals with similar statistics (as the C/A codes are)? If this is not the cause, what do you think it is?
Thanks in advance.


texts: QPSK demodulation problem due to apparent transmitter clock frequency drift

texts:
I am using FPGA to generate 10Mbaud QPSK signal in two channels corresponding to two quadratures. Upsampling to 1GSample/s (using raised cos filter for interpolation) and complex upconversion (f_up=250MHz) is done in FPGA. This signal is sent to a DAC and detected back-to-back using Tektronix rsa607a spectrum analyzer. I sent 200k sampes and detected them using the analyzer.

This is a screenshot form the measurement.  On the right side is unzoomed waveform in time where red dots represent sampling points in demodulation. On the left is constellation for the whole stream. As you can see the analyzer synchronizes well with the middle of the captured signal, however it goes off sync with time. If I measure shorter intervals EVM gets smaller. On the left side you can see measured frequency error which is the difference between transmitter frequency and analyzer's internal clock. This error is obviously different with time which tells me the FPGA clock might be unstable. However, I did other measurements when offloading a sine tone where I don't see these frequency drifts.
Does anybody have intuition what can be the cause of this problem?

Cheers,
D


texts:
As a part of my thesis work -in which I am using both acoustic and visual information- I am looking forward to implement a real time sound source localization sub-system. From the geometric principles, I know that at least 4 non-coplanar microphones are needed for 3D localization.

However, as I understand, using only a laptop to interface 4 microphones is impossible. I read that one sound card can only support 2. I also search through 4 microphone usb interfaces and found out quiet expensive solutions. I have also encountered to PS Eye, I'm not sure its sensory information is accessible through MATLAB or CLAM/C++ though.

As I only need an approximate information, I then decided to use only the azimuth angle information (planar confidence) of the source. For that, I know 2 microphones are enough. But I read about another theoretical problems for this case, like the need of calibration, harder sound seperation etc.

I have to ask is there a cheap interfacing solution, using 3 or 4 USB microphones? If not, would it be feasible to build up something from scratch and use amplified electret mics? I am really confused about hardware selection. Neither my advisor, nor the academic papers have too much practical information.

Thanks in advance.


texts:
I understand why most of the classic textbooks only consider Additive White Gaussian Noise channels in digital data transmission. Unfortunately in my case noise is not Gaussian and it might even be colored (non-white).

I spent a couple of hours searching around but only found some special cases that were eventually approximated by Gaussian models. I appreciate if you could point me to some resources that non-Gaussian noise is considered, in concise and direct manners.

Some notes: I'm interested in evaluating the performance of a receiver in the presence of Gaussian and other non-Gaussian noise sources. For example, I'd like to see the effect of the ADC quantization noise passed through some filter, like an equalizer (and there is not much support to use Central Limit Theorem to make Gaussian assumption), or some residual jitter from my clock recovery loop that again doesn't fit nicely to a Gaussian model.


texts: How to Simulate AWGN (Additive White Gaussian Noise) in Communication Systems for Specific Bandwidth

texts:
I am trying to generate a AWGN waveform to add it to the signal of my simulated communication system. The operating bandwidth of the communication system is about B=3GHz and suppose that T=300K (my signal is comprised of very short pulses, I am simulating an IR-UWB system).

What I have done so far is:

N = k*T*B
sigma = sqrt(N./2)
noise = sigma.*randn(1,s_length)

where s_length is the number of samples for my useful signal and noise the awgn waveform.

Is this procedure correct? I have browsed through different questions but was unable to clarify it.


texts:
I have implemented 256-QAM based OFDM design on FPGA. The data flow is as follows -
QAM -> Hermitian Symmetry -> IFFT -> Pilot Insertion -> Cyclic Prefix
I have done Hermitian Symmetry before the IFFT block, so that the imaginary part becomes zero. Now, the O/P of Cyclic Prefix block is transmitted to DAC with some additional known pattern padded at the beginning of each frame for synchronization at the receiver end. The output of DAC is directly connected to the ADC input via Co-axial Cable, having amplifiers and LPF in between. The specifications of the system are as follows -
OFDM Signal Frequency - 40 MHz
DAC Resolution - 12 bits
DAC Sampling Rate - 200 MSPS
ADC Resolution - 12 bits
ADC Sampling Rate - 200 MSPS
Problem Statement - I am unable to recover the original values that were transmitted to DAC. The issue is that the deviation of ADC Output is non linear as compared to the DAC Input values.
How to recover the original values? I have implemented the Linear Interpolation based Channel Estimation also, but it is not able to compensate the error cause by hardware noise. Until n unless the correct values are fed to the FFT block, the de-mapped data output of De-QAM won't be correct.
Need help in recovery of the OFDM signal.
P.S. - PFA the waveform snapshot of the Transmitted/Received OFDM Signal. Blue waveform represents Transmitted Signal; Orange waveform represents Received Signal.

Thanks


texts: Proof of Paley-Wiener criterion for causality

texts:
The Paley-Wiener criterion for causality is that $\displaystyle\int_{\mathbb{R}}\frac{A(\omega)}{1 + \omega^2}\mathrm{d}\omega$ exists and is finite, where $A(\omega) = \left|\mathcal{F}[f]\right|$ is the modulus of the Fourier transform of a function $f: \mathbb{R}\to\mathbb{R}$

How is this derived from the usual statement of the Paley-Wiener theorem?

I have also seen that this proves that a causal filter can't provide total rejection for a band of frequencies, but this appears to me to follow directly from the usual statement that the Fourier transform of a square-integrable distribution vanishing for negative reals is holomorphic and the identity theorem for holomorphic functions.


texts: get convex hull of isolated shapes with openCV/C++

texts:
I have segmented CT images, something like the following one:

I would like to get a convex hull that contains all the white regions.
Something like this circle, but more fitting, and of course not necessarily a cicle.

I've already tried this example https://stackoverflow.com/questions/11813352/finding-the-convex-hull-of-an-object-in-opencv but I've only got this output:


texts: Unsupervised Way to Choose K in K-Means Algorithm

texts: There are couple of methods listed here about finding the value of k in k means algorithm. How much sound is each of these procedure when used to find out the value of k in unsupervised manner.


texts: simulator for dsp kit tms320c6713?without using actual hardware kit tms320c6713?

texts: As you know, most universities across the globe are closed due to pandemic  and people are unable to acess research facilites ,so in such a scenario one wants to improve his dsp hardware programming skills, how can he make up for dsp kit tms320c6713? I mean, is there any software proteus/matlab or Labview or anyother software that enables to simulate tms320c6713 without using actual hardware tms320c6713


texts:
In the book Modern Digital and Analog Communication Systems by B. P Lathi, there is a section for the estimation of WBFM bandwidth. First, the staircase approximation of $m(t)$ is constructed by the pulses of width $1/2B$. $$\hat{m}(t) = \sum_{k=-\infty}^{+\infty}m(t_k)\Pi(2B(t-t_k))$$
It's stated that FM signal corresponding to each of this cells is a sinusoid with frequency $\omega_c + k_fm(t_k)$ with duration $1/2B$. Why the corresponding FM signal is of finite duration? After that it's stated that "The FM signal for $\hat{m}(t)$ consists of a sequence of such constant frequency sinusoidal pulses of duration $T = 1/2B$ corresponding to various cells of $\hat{m}(t)$". It seems that constructing the FM signal assumed to be linear here. Or maybe I'm missing something?
Source: http://www.cs.csub.edu/~vvakilian/CourseECE423/LectureNotes/Lecture9.pdf
Edit: Actually, there are two main questions here:
1. Why the FM signal corresponding to $m(t_k)\Pi(2B(t-t_k))$ is a sinusoid with frequency $\omega_c + k_fm(t_k)$ with duration $1/2B$? (I think it should be infinite duration.)
2. Why the FM signal corresponding to $\sum_{k=-\infty}^{+\infty}m(t_k)\Pi(2B(t-t_k))$ is sum of the sinusoids with duration $1/2B$ and frequencies $\omega_c + k_fm(t_k)$? (Assuming that FM signal corresponding to each cell is is a sinusoid with duration $1/2B$, it's not true that FM signal corresponding to sum of these cells is equal to sum of the FM signals corresponding to each cell.)


texts:
I have a signal h that I want to delay for a time t.
I know that I could use two possible approaches:

$h_{\text{delay}}[k] = \delta[k-\tau]\star h[k]$
$h_{\text{delay}}[k] = \text{IFFT}\left(\lvert H(\omega)\rvert e^{j\phi(\omega)}\cdot e^{-j\omega \tau}\right)$

Which theoretically should be completely equivalent (as far as I know).
From the performance and precision point of view, which of the 2 is better? Is there any trade-off in numerical precision or speed?


texts: ECG SIGNAL ANALYSIS USING ADAPTIVE THRESHOLD

texts: I am trying detect the R-peaks in ECG signal using adaptive thresholding method in MATLAB. I have done till preprocessing of the signal and find the local maxima of the signal but I am unable to understand how to set the adaptive threshold from those points, anyone has an idea about how to do it?


texts:
After investigating signal analysis methods based on empirical mode decomposition (EMD), I found that recent developments are mostly related to the Hilbert Huang transform (HHT) and the Local Mean Decomposition (LMD) method.

I have been reading a few articles on the subject and I would like to have your opinion regarding HHT.

EMD leads to intrinsic mode functions that are particularly well suited for Hilbert transform. HHT seems to be really widely used for many distinct types of industrial or academic applications.

Am I right to think that HHT is essentially EMD + Hilbert transform? Do you consider the specificity of the HHT lies in the sifting process based on spline interpolation (to the contrary of the LMD using a moving average algorithm for instance)?


texts:
I've created both single tone signal and a LFM chirp signal with same pulse duration and cross correlated both of them with their own(matched filter). Peak points of both results were same as expected which is the total energy of input pulse. I've read about SNR gain with signal compression in radar equations and chirp is an example of compressed signal. Where does this SNR gain come from? Does it simply mean matched filter? If so pulse compression is no different than a single tone as far as I understand(range resolution is the only difference).

Edit:
Equation 7.4 states that output peak SNR to input SNR is proportional to both bandwidth and pulse length. It looks at the output to input SNR ratio which is affected by the input noise because of bandwidth used at the input. Does this mean a single tone signal has a better input SNR because of less bandwidth? Then what is the paper trying to say by comparing output SNR to input SNR if peak point is the total energy of my pulses on both compressed and single tone after a matched filter?

Edit2: After some thinking I've concluded it as follows: No matter what the signal modulation or bandwidth, peak SNR at matched filter output is same for input pulses of same energy. However, input noise power is proportional to the bandwidth of signal which gets filtered on matched filter and that's what we call matched filter gain. Correct me if I'm wrong.


texts: Are there any recursive online max/min filters for time-series

texts: Are there any online recursive filters that can approximate local, time-varying minimum and maximum values of a time series?


texts: How do I find a system's impulse response from its state-space repersentation using the state transition matrix?

texts:
Suppose we have a linear represented in the standard state space notation:

$$ \dot{x}(t)=Ax(t)+Bu(t)$$
$$y(t) = Cx(t) + Du(t)$$

In order to get its impulse response, it is possible to take its Laplace transform to get

$$sX=AX+BU$$
$$Y=CX+DU$$

and then solve for the transfer function which is

$$\frac{Y}{U}=C(sI-A)^{-1}B+D$$

Similarly, for a discrete system, the $\mathcal{Z}$-transform of
$$ x[n+1]=Ax[n]+Bu[n]$$
$$y[n] = Cx[n] + Du[n]$$

is

$$\frac{Y}{U}=C(zI-A)^{-1}B+D$$

This process seems a bit long and I remember that there is a way to find the impulse response using the state transition matrix which is the solution for $x$ of the first equations of each pair. Does anyone know how to do this?


texts: Unit impulse response of a cascade interconnection of three discrete-time systems

texts:
I am nearly at the end of finishing a problem in my textbook but I couldn't understand something in the answer;

I did everything to the point I found the overall response of the system in terms of $h_{1}[n]$ which is $$h[n]=h_{1}[n]+2h_{1}[n-1]+h_{1}[n-2]$$ Now while comparing each $n$ value with the given graph for overall response, the answer is saying that for $n=0$;  $$h[0]=h_{1}[0]$$ but I didn't get why it isn't like this: $$h[0]=h_{1}[0]+2h_{1}[-1]+h_{1}[-2]$$ How do we know $h_{1}[n]$ is zero for $n$ values smaller than zero?


texts:
Let's say that we have two sequences, input sequence $x(n) = [0121]$ and impulse response of a given system $h(n) = [0, 1, -1, 1]$. I need to find response of this system to given input sequence. After that, i need to calculate linear convolution of given sequences.

If we denote the response as $y(n)$ we have $y(n)=h(n)*x(n)$, which means that, due to convolution theorem, in frequency domain we have $Y(k)=H(k)X(k)$. From this, we can find $y(n) = IDFT(Y(k))$.

Considering the fact that i need to find convolution of given sequences, meaning $y(n)=h(n)*x(n)$, that would mean that it should yield same result as when i was doing this by using DFT. However, my final results don't match at all.

$x(n) = [0121] \Rightarrow X(k) = [4, -2, 0, -2] \\ h(n) = [0, 1, -1, 1] \Rightarrow H(k) = [1, 1, -3, 1] \\ Y(k)=X(k)H(k) = [4, -2, 0, -2] \Rightarrow y(n)=[0, 1, 2, 1]$

On the other hand, convolution of given sequences gives the following result:

$y(n) = h(n)*x(n) = [0, 0, 1, 1, 0, 1, 1]$

Not only that result is completely wrong, but also dimensions of vectors i got as a result are not the same. What am i doing wrong? Any help appreciated!


texts: Why DTFT coefficients are periodic and why continuous Fourier transform coefficients are not periodic?

texts: As I understand, when the input signal is discrete in time and we want to find the coefficients of Fourier transform, DTFT is used and the coefficients in frequency domain are periodic, but I can't understand why? If DTFT coefficients are periodic, why the coefficients of a continuous signal in frequency domain are not periodic?


texts: Is Image Noise Removal Image Enhancement or Image Restoration Operation?

texts:
I am working in digital image processing field, recently I am studying about Image Noise, I just want to know that whether this noise removal is part of image enhancement or image restoration.

I have read some papers on it, at some points authors are saying that image enhancement from a given noisy image, and on the other hand in some papers authors are saying that noise removal is image restoration process.

Please explain what is correct?.


texts: how do I derive the system equation for a simple delay with feedback?

texts:
I am a software engineer, and just learning digital signal processing formally, though I've hacked around before a fair amount.

I'm implementing a delay audio VST and I'm trying to wrap my head around how one creates equations of the output given a system diagram. I'm not used to circuit diagrams (they're weirdly parallel to me), so I'm trying to figure out how one goes from this system diagram:

to an equation. As best I can tell, the equation for this is:

y[t] = x[t] + y[t - T] * gain

where x[t] is the input on the far left, y[t] is the output on the far right, and the delay time in samples is T.

Is this correct?

I don't fully see how the system diagram implies the second term is y[t - T] rather than x[t - T], but intuitively, I know it must be, because an audio delay unit can't really incorporate feedback if the original signal delay only happens once (feedback must enable the ever-decaying ringing with time). Otherwise we'd get just a shifted version of the original, added to the original. This would be boring.


texts: Train Pulse Recognition in Medium Voltage Networks

texts:
I just started a project in the Smart Grids area. It involves the recognition and quick extinction of earth faults in a medium voltage system. The system works automatically and pretty much gets rid of earth faults very quickly.

My task is to make and algorithm to also know where the earth fault happened along the network. For that, the system is programmed to send a rectangular pulse train in the line, which modulates the neutral current (which is close to zero much of the time). Once the current is modulated with the pulse train, the work I have to do is to recognize the pulse train by analyzing the neutral current. Once the pulse train is recognized at a given substation I know where the fault happened.

The pulse train that is sent is always the same and has known amplitude, frequency and number of pulses (6 amperes, t_on=t_off=1s and 5 pulses, for example) and cannot be changed for this project. Right now the system has an algorithm running to locate the fault by doing an FFT on each pulse it encounters and measures the harmonics energy in order to see if it is in fact a rectangular pulse or not (since the FFT of a rectangular pulse has a know spectrum). Once the five pulses are positively recognized one after the other, the Earth Fault Detected flag is raised.

The problem is that outside pulses can also come on the line randomly, from say a wash machine being turned on and off or a blender, for example. This leads to false earth faults being detected.

My initial solution and recommendation is to use the cross-correlation and auto-correlation function on the whole pulse train in order to check periodicity of the signal. Once that's done I can know if the signal does in fact have this particular pulse train or not. I did a small Matlab simulation to see the ACF of my pulse train and of a random set of pulses to further cement my hypothesis:

Cross Correlation:

Autocorrelation:

Does this though process make sense? Does anyone have any other recommendations to accomplish this?


texts: pairwise error probability and union bound of space time code

texts:
i have a question that might seem silly.
For a given space time code, say a constellation of 16 unitary matrices. When using the union bound to derive an upper bound on the average error probability (using the equation in the image below)

is it logical to have the Bit Error Rate result higher than 1 for a given SNR? even if what i am deriving is an upper bound ?

Thank you


texts: How to replace Hadamard Product and Column wise addition with new Matrix operation?


texts: How to obtain a meaningful group delay graph using matlab

texts:
Problem:
I'm trying to analyze the behavior of an FIR filter with the following impulse response/kernel:

Using Matlab's function grpdelay(myKernel,length(myKernel)), I obtained the following figure:

Based on my limited knowledge the graph is not correct since it contains negative values. Is my understanding correct? If it is, how could I go about obtaining a more accurate group delay graph?

Code:

    onePulseRun20 = dlmread("C:\Users\agusfrpa\Pictures\singlePulseRun20.txt");
twoPulseRun21 = dlmread("C:\Users\agusfrpa\Pictures\twoPulseRun21.txt");
%twoPulseRun22 = dlmread("C:\Users\agusfrpa\Pictures\twoPulseRun22.txt");

%y is averaged one pulse data
y = mean(onePulseRun20,1);
pulse = twoPulseRun21(1,:);

myFilter = conj(fft(y))./(abs(fft(y)).^2+.3);
myKernel = real(ifft(myFilter));
myKernel = myKernel(4700:5250);

grpdelay(myKernel, length(myKernel));

I will post the link for data in the comments.


texts: How to find roll-off from its transfer function equation?

texts:
I want to find roll-off (unit dB/decade) of highpass filter $H(z) = 1 - a \cdot z^{-1}$.

I derived magnitude response to be

$1+a^2 - 2a\cdot \cos(\omega)$

To find out roll-off factor from above equation, how should I proceed?

I thought of doing

$1+a^2 - 2a\cdot \cos(10\omega) \over {1+a^2 - 2a\cdot \cos(\omega)}$

But $\omega$ is a variable, how do I simplify this?


texts: Finding time of arrival of signal using GNU radio

texts:
I am sending a repeating BPSK modulated PN sequence signal from a USRP x310 transmitter. I am receiving the sequence with another USRP x310. Both are attached to an octoclock. Is there a way to get the global time when the receiver starts writing to the file sink?
Here is a look at the receiver block diagram:

Currently, I can use the file metadata sink to get the time it starts writing (starting at 0 sec) and subsequent writes after a number of samples. But it does not give the global time. I eventually want to use this global time measurement to calculate the time difference of arrival at multiple receivers.


texts: Synthetic Aperture Radar (SAR). How to do Range Compression?

texts:
I am trying to perform range compression to a raw SAR image [fast times, sensor position], Dr in code. I have the transmitted signal (chirp), g in code, and literature says range compression is obtained by convoluting the time reverted conjugate of g along the rows. I use the following MATLAB code.

    h_compressed = conv2(conj(fliplr(g), 1, Dr, 'same');

Where 3003 is the number of columns of the RAW image. Nevertheless my output is the same as my input i.e. no range compression occurs. Any ideas on why?

After range compression I perform 1D matched filtering.

    for k = 1:Nr
    R = sqrt(x_filter.^2 + SAR_data.r_ax(k).^2);
    filter = exp(-1i*4*pi/lambda*R);
    Dfr(:,k) = conv2(Drc(:,k),filter','same');
    end

And obtain the following:


texts: Bandpass Experiment: Weird graph! amplitude mismatch?

texts:
Why is this generic BP filter seem to magnifying amplitudes when I run it at lower sample rates? The bass range seems to 32 times louder at 1275 Hz (44.1 kHz / 32) than at 44.1kHz?

Here is a graph of a birdsong which has been passed through 128 BandPass filters on MIDI notes 0 to 128.

To make the graph draw faster, I have attempted an experiment: skip 32 audio samples for frequencies lower than 150HZ.

The filter-bank draws 15 times faster, but the results are wrong: the frequencies lower than 150Hz seem to have been multiplied by an amplitude of 32? What is happening?

correct graph:

The same downsampling method doesn't amplify a sine sweep, so I think it will not be possible to use this computation optimizing trick for a high quality spectrogram?


texts:
I have darkcurrentrate (Dark noise) data collected at various ISO speeds such as ISO 55, ISO 99, ISO 198 etc. and the files are in .dng format.
Can someone suggest how to calculate the dark current rate, which can be derived by checking the slope of the dark image with different exposure times using MATLAB?


texts: How composite signal in CDMA is normalise in general?

texts:
I had CDMA composite signal with 3 users. I assume that the signals are linearly sum together so for my maximum and minimum amplitude of composite signal are +3 and -3 respectively.

How can i normalize my signal, generally its by number of users at receiver side. But for transmission, don't i need signal with +1 and -1 ? Am i missing something ? What would be the correct approach ?


texts: Compression algorithms specific to complex signals

texts:
I am looking for (lossy or lossless) compression algorithms dedicated to complex signals. The latter could be composite data (like the left and right for stereo audio), a Fourier transformation or an intermediate step of a complex processing, an Hilbert pair of generic  signals, or complex measurements like some NMR (Nuclear magnetic resonance ) or mass spectrometry data.

A little background: I have been working on real signal and image coding, and mildly on video and geological mesh compression. I have practiced adaptive lossless and lossy compression (RLS, LPC, ADPCM), Huffman and arithmetic coding, transforms (DCT, wavelets), etc. Using standard compression tools separately on each of the real and the imaginary parts is not the main aspect of the question. I am interested in compression-related ideas specific to complex data, for instance:

sampling of the complex plane,
joint quantization of module and phase,
optimal quantization of 2D complex values: Loyd-Max is "optimal" for 1D data. I remember that 2D optimal quantization is generally more complicated. Are there 2D binnings dedicated to  complex, or for instance Gauss integers?
entropy coding methods, arranging along complex "features" (angle, magnitude),
entropy calculation for complex data,
use of specific statistical tools (e.g. from Statistical Signal Processing of Complex-Valued Data: The Theory of Improper and Noncircular Signals)
integer transformations for Gauss integers.

I have not found much tools. My scarce references are:

An experimental audio codec based on warped linear prediction of complex valued signals, 1997
Compression of Complex-Valued SAR Images, 1999  (Download)
Lossless Signal Processing with Complex Mersenne Transforms, 2003

I would be interested in more authoritative references. A compression file format dedicated to complex data would be a plus.


texts: Sampling: How should I sketch a signal spectrum?

texts:
As revision, I am going through Signals and Systems, and in Section 7.1the authors are sketching representative spectra for x(t) and xp(t)as follows

However when I search for "How to draw a signal spectrum" I get results (here and here), where signal spectra are shown like this:

For example,how should I go about sketching the the spectrum of signal x(t), where x(t) = 4cos(8000πt)?

Here are my two answers:

Which one is the correct approach? Any suggestions would be appreciated.


texts:
I'm studying on designing 9-tap zero forcing equalizer to remove ISI in the channel. I assumed that the channel response is $h[n] = [0.41, 0.815, 0.41]$ and I tried to calculate bit error rate. When I choose $h[n]=[0.21, 0.815, 0.21]$ for example, I get good results. BER goes to 1e-4 for high SNR values. But the code below gives bit error rate of 0.5 when the channel response is $[0.41, 0.815, 0.41]$. Can anybody help me about this problem? I did something wrong but I couldn't find it.
N  = 10^4; % number of bits or symbols
Eb_N0_dB = [0:2:10]; % multiple Eb/N0 values

for ii = 1:length(Eb_N0_dB)

    % Transmitter
    ip = rand(1,N)>0.5; % generating 0,1 with equal probability
    s = 2*ip-1; % BPSK modulation 0 -> -1; 1 -> 1

    % Channel model, multipath channel
    ht = [0.41 0.815 0.41]

    chanOut = conv(s,ht);
    n = 1/sqrt(2)*[randn(1,N+length(ht)-1) + j*randn(1,N+length(ht)-1)]; % white gaussian noise, 0dB variance

    % Noise addition
    y = chanOut + 10^(-Eb_N0_dB(ii)/20)*n; % additive white gaussian noise

    kk=4; % 9-tap filter
    L  = length(ht);
    hM = toeplitz([ht([2:end]) zeros(1,2*kk+1-L+1)], [ ht([2:-1:1]) zeros(1,2*kk+1-L+1) ]);
    d  = zeros(1,2*kk+1);
    d(kk+1) = 1;
    c  = [inv(hM)*d.'].';

    % matched filter
    yFilt = conv(y,c);
    yFilt = yFilt(kk+2:end);
    yFilt = conv(yFilt,ones(1,1)); % convolution
    ySamp = yFilt(1:1:N);  % sampling at time T

    % receiver - hard decision decoding
    ipHat = real(ySamp)>0;

    %     % counting the errors
    nErr(1,ii) = size(find([ip-ipHat]),2);
end

simBer = nErr/N; % simulated ber


texts: Standard representation of FFT of an Image?

texts:
Statement 1 : When we take DFT of an image matrix, we get DFT co-efficient matrix. What does  the DFT Co-efficients tells about the image? (For example, take the 1st co-eff in the DFT matrix, what does it tells us ?)

Statement 2: In a DFT matrix, we have low frequency components on the top left corner of the matrix and the frequency will be getting increase as we move down and we have high frequency in the bottom right corner. Is this statement correct ? Correct me if I am wrong.

I couldn't able to understand the following image (which contradicts my statement 2 (I think so )) , can anyone explain this pls?


texts: Data fusion using 2d discrete wavelet transform (DWT)

texts:
I am working on a project that employs a linear array sensor that provides data from the same object at two different energies.  Collected in time, I end up with two images (16-bit sensor values, MxM image), call one HIGH the other LOW.  My engineering background is Material Science and metallurgy, so I am coming up to speed on DSP, image analysis, and all of the bits that bring all of that together.

I have been researching "image fusion" and have been working from Gonzalez's fabulous Digital Image Processing (3rd ed.) text.  I have developed an understanding of DWT and have successfully written code to replicate the 3-scale DWT transforms in the text (major accomplishment for me).

Now, I am trying to understand how to "fuse" the DWT coefficients from the LOW and HIGH images.  From the IEEE Xplore database, I pulled several papers on the topic; one of the more useful is by Zhiyu Chen, et.al. "A Combinational Approach to the Fusion, De-noising and Enhancement of Dual-Energy X-Ray Luggage Images".

Succinctly, they recommend "averaging the corresponding approximation coefficients of L and H" (sigma)... and "summing the corresponding detail coefficients of L and H" (psi).  But, that is all of what is said on the matter.

I have several other sources, but they must presume some obvious bit of understanding that I am lacking!

My question is, in order to successfully fuse the HIGH and LOW images: what does this mean?  In the case of a 3-scale DWT, do I average or sum the 3-rd scale approximation/detail coefficients?  Do I somehow average across each scale, despite the fact that the size of the coefficient matrices are different by powers of 2?


texts:
I have been using Garmin GPS device for years. Just curious which algorithm is used in the device:

Navigation equations that can calculate both GPS device's position and its clock bias. In this way, the time display should be accurate with the clock in the satellite. However, the computational complexity could be very high. See https://en.wikipedia.org/wiki/Global_Positioning_System.
Use TDOA method similar to what is explained at https://en.wikipedia.org/wiki/Multilateration. The computational complexity could be much lower, but you will not be able to obtain the receiver's time bias.

Note: This question is based on the assumption that the satellites' positions and TOFs have already been extracted from the received satellite signals.


texts:
I am generating an AM modulated sinusoidal wave. The carrier frequency is set at $1000 \ \rm Hz$; the modulation frequency is set to $40 \ \rm Hz$ and $100 \%$ of the amplitude is modulated.
I generated 2 signals with 2 slightly different equations.

Signal 1 follows:
$$\sin(2\pi f_m t) \cdot \sin(2\pi f_c t)$$
while signal 2 follows:
$$\big(1-\sin(2\pi f_m t)\big) \cdot \sin(2\pi f_c t)$$

I don't understand why the first signal doesn't have a $1 \ \rm kHz$ component on the FFT and why both signals do not have the same period.
Both $\sin(x)$ and $1-\sin(x)$ have the same period, and I can see that signal 1 has the amplitude sinus and the carrier sinus in phase, but I can't get my head around those plots and interprete them correctly. Thanks for all the information and explanation you can provide :)

Code snippet to produce the 2 signals:
#!/usr/bin/env python3
# -*- coding: utf-8 -*-

from matplotlib import pyplot as plt
import numpy as np
import os

#%% Signal 1
fs = 44100        # sampling rate, Hz, must be integer
duration = 1.0    # in seconds, may be float

# Define the time series
t = np.linspace(0, duration, int(duration*fs), endpoint=False) # time variable

# AM - Amplitude Modulation
fm = 40 # Modulation frequency
amplitude = np.sin(2*np.pi*fm*t)

# Output signal
fc = 1000 # Carrier frequency
signal1 = amplitude * np.sin(2*np.pi*fc*t).astype(np.float32)

# Apply FFT
fft_freq1 = np.fft.rfftfreq(signal1.shape[0], 1.0/44100)
fft1 = np.abs(np.fft.rfft(signal1))

#%% Signal 2
fs = 44100        # sampling rate, Hz, must be integer
duration = 1.0    # in seconds, may be float

# Define the time series
t = np.linspace(0, duration, int(duration*fs), endpoint=False) # time variable

# AM - Amplitude Modulation
fm = 40 # Modulation frequency
amplitude = np.sin(2*np.pi*fm*t)

# Output signal
fc = 1000 # Carrier frequency
signal2 = (1-amplitude) * np.sin(2*np.pi*fc*t).astype(np.float32)

# Apply FFT
fft_freq2 = np.fft.rfftfreq(signal2.shape[0], 1.0/44100)
fft2 = np.abs(np.fft.rfft(signal2))

#%% Plot
f, ax = plt.subplots(2, 3, sharex=False)
ax[0, 0].plot(t[:4411], signal1[:4411])
ax[0, 0].set_title('Signal 1')
ax[1, 0].plot(t[:4411], signal2[:4411])
ax[1, 0].set_title('Signal 2')

ax[0, 2].plot(fft_freq1[900:1101], fft1[900:1101])
ax[0, 2].set_title('Signal 1 FFT')
ax[1, 2].plot(fft_freq2[900:1101], fft2[900:1101])
ax[1, 2].set_title('Signal 2 FFT')

ax[0, 1].plot(t[:4411], amplitude[:4411])
ax[0, 1].set_title('Signal 1 AM')
ax[1, 1].plot(t[:4411], (1-amplitude)[:4411])
ax[1, 1].set_title('Signal 2 AM')


texts: Lets say I have a small size vector x=[a b c d]. Now I stretch this vector 3 times and I got x3=[a a a b b b c c c d d d]. What would be the relation between fft(x) and fft(x3)?


texts: Reconstruction of a signal from non-uniform samples

texts:
So I'm having some troubles with a signal that's been sampled but where the samples don't have consistent intervals. I've been looking into the NFFT (in python) as a starting point but I'm a complete noob with signal processing so I'm not really getting anything out of it.

So I have one array with the sample values and another with the corresponding timestamps. But I really have no idea how to get the spectrum out of it. I've done the first part as described here.

Sorry about the sparse question.


texts: How to estimate the autocorrelation from nonuniformly spaced data

texts:
Assume a continues-time random process $X(t)$ sampled nonuniformely in time to acquire discrete signal $x[n]$. The sampling times are known but the autocorrelation is not. Is there an accurate approach to estimate the autocorrelation function? This problem is challenging since at some lags there is no sample to average.

Any suggestion even in a special case, e.g. a simple autocorrelation function with an exponential decay (continuous-time AR process of order one) would be appreciated.

Edit

I have already read some papers, so I would like either an analytical answer or a peace of algorithm/code to address this problem for instance for an AR(1) process assuming non-uniformly spaced samples(e.g. in uniformly distributed random times).


texts: transient response of simple digitalized RC low-pass filter

texts:
In http://en.wikipedia.org/wiki/Bilinear_transform#Example, digital version of simple RC low-pass filter is presented:

$$\frac{1 + z^{-1}}{(1 + 2RC / T) + (1 - 2RC / T) z^{-1}}$$

where $T$ is sampling interval. Let me reorganize by $A = 2RC/T$, and get the following:

$$y(n) = \frac{x(n) + x(n-1)+(A-1)y(n-1)}{A+1}$$

But I cannot see how I would get transient and steady response back from this difference equation. Can anyone help here?


texts: How does the error in azimuth angle vary as one moves away from the boresight direction of a radar

texts: My question is basically same as the title. I am trying to find a mathematical model which provides the error, $d\theta$, in azimuth angle($\theta$) as one moves away from the boresight direction of radar. It is also known that the error in the radar data is also a function of distance ($r$) from the radar. So, in order to incorporate the distance from the radar, can anyone tell me how model $d\theta$ as a function of $\theta$ and $r$?


texts:
This is now a second time I am attempting to ask this very important but simple question here. What I want to know is can you do deconvolution by convolving a signal. It is often stated that, for example by cutting and boosting the same frequency on an equalizer the result is the original signal. Is that the case? Can convolution be removed by convolving? That would certainly brake the identity that the convolved signal must be n+m-1 in length. I tried this with an equalizer and the results seem to be close to perfect, just some quantization distortion.

If the above can be done, I can not arrive at this conclusion. For example If I convolve (1) with (1, 1) the result is (1, 1). I can not find any impulse response that if convolved with (1, 1) would result in (1).

So what is the truth here?


texts: In the Sampling Theorem, why are the image frequencies at n*fc not a problem?

texts:
In this example I have been working through, we first look at the situation when fs > 2fc and then the situation when it isn't:

In the example, the frequency responses of a sampled sinusoid at fc = 5 kHz are shown. In the first instance, the sampling rate is fs = 16 kHz, and in the second it is fs = 8 kHz.

I understand the concept of aliasing and I can see how it has occurred in the second example - an alias has crossed over fc. What I don't quite get is why the situation in the first instance is acceptable - there are still aliases in the signal. Is it simply that we can window (low-pass filter) the frequency response to fs/2 as we know this is the highest frequency that can be resolved?


texts: Why is a precoder necessary for DQPSK and what does it accomplish?

texts:
I've implemented a soft-decoder for DQPSK using the wonderful answers I received here:

How to soft decode DQPSK?

To get the soft-decoder working properly I needed to precode the data I was sending out.  I implemented the precoder mentioned in this paper:

$I_k=\overline{u_k \oplus v_k}*(u_k \oplus I_{k-1})+(u_k \oplus v_k)*(v_k \oplus Q_{k-1})$
$Q_k=\overline{u_k \oplus v_k}*(v_k \oplus Q_{k-1})+(u_k \oplus v_k)*(u_k \oplus I_{k-1})$

I'd like to know why this precoder is necessary -- what does that complicated expression of XORs actually accomplish?

Here's a table showing what the equation yields.  If "to_encode" is 00, the to_send symbol is the same as the previous ("prev") symbol.  If the "to_encode" is 11, the to_send symbol is the previous symbol xor 11.  What is the meaning in other cases?

to_encode prev  to send
[ 0 0 ] [ 0 0 ] [ 0 0 ]
[ 0 1 ] [ 0 0 ] [ 1 0 ]
[ 1 0 ] [ 0 0 ] [ 0 1 ]
[ 1 1 ] [ 0 0 ] [ 1 1 ]
[ 0 0 ] [ 0 1 ] [ 0 1 ]
[ 0 1 ] [ 0 1 ] [ 0 0 ]
[ 1 0 ] [ 0 1 ] [ 1 1 ]
[ 1 1 ] [ 0 1 ] [ 1 0 ]
[ 0 0 ] [ 1 0 ] [ 1 0 ]
[ 0 1 ] [ 1 0 ] [ 1 1 ]
[ 1 0 ] [ 1 0 ] [ 0 0 ]
[ 1 1 ] [ 1 0 ] [ 0 1 ]
[ 0 0 ] [ 1 1 ] [ 1 1 ]
[ 0 1 ] [ 1 1 ] [ 0 1 ]
[ 1 0 ] [ 1 1 ] [ 1 0 ]
[ 1 1 ] [ 1 1 ] [ 0 0 ]


texts: Splitting of signal into different components for calculation of PSD

texts:
In BCI competition 3 dataset 5, It has been stated that we have to take 8-30 Hz frequency and break it with resolution of 2 Hz then it will result in 12 components but I am finding out only 11 components.Why??
The description can be found out here...

http://www.bbci.de/competition/iii/desc_V.html


texts: Frequency at which the maximum amplitude occurs of a signal

texts: How can i  find the frequency corresponding to maximum amplitude of each frame after uniform segmentation.


texts: How to build a software audio mixer - low level approach

texts:
my questions might come across different SE forums. But here it is.

A little context : I'm currently building a midi sequencer, the architecture is this way (for now) I've got a ARM v7 + a PIC, the ARM v7 runs a linux distro and will do the Human Interface (getting interactions on buttons and hardware interfaces + react to them in order to fill my loops) and the PIC mcu, will only do the proper sequencing, ie : it will only write in rhythm what the arm will tell him. For now, the pic is only doing the real time interface with MIDI out using a crystal, at some point it might also handle some other hardware feature such as getting the inputs and outputs from the buttons and rotary but I will decide that later. And the linux board might only do audio processing (sampler part) + UI.

Now what I am interested in is the following :
I'd like to add into my system a sampler, the idea is to have up to 16 samples and to play them into an audio output. I've already used some higher level approach to do that using external libraries as FMOD or even Alsa or Gstreamer, but what I want to do is to have full flexibility and possibility of optimization on my mixer. For instance later part will include effects and tuning, and bit crushing, and this kind of stuff. In the first version I'd like to add the amplitude of the signal, ie a value between 0 and 1 that would set the amplitude of the sample regarding the state of a rotary encoder.

So here is my question, let's say I have 16 wav file, 48kHz, 16bits. How can I do the mixing of these 16 samples ?

If I call k[i] the amplitude setting, and samplei the pcm data of the sample at the current moment then is the following correct :

mixer_out(t) = (k[0]*sample [0] (t)+...+k[15]*sample [15] (t)) / (k[0]+...+k[15])

I've seen every thing and it's contrary when it come to sample mixing. First can I simply add my pcm data in order to get the mix of the sample ? Then (if we don't take into account the clipping) how shall I afterward put my sample mixed into the mixer_out ?
I think my naive idea will be far from a correct way of doing the mixing, but what I need is a very low level (at least at pcm level) approach of mixing.

So let's take N pcm samples (all at the same format), and 1 pcm audio out called mixer_out, and let's mix all the sample together adding a constant k[i] (i from 0 to N) for all the samples (k is a float between 0 and 1).
I believe that it could be possible to take k from 0 to 1/N in order to simplify if that helps. Even though I do understand that this won't be optimal when I will have several sample with k=0.

Also take into consideration that my system will be run into an embedded processor in "real time" (or do be more precise, the mixer out will have to be calculated at least as fast as or faster than the audio is played) so it will not have unlimited processing power, but I'm ready to do compromise with the number of output in order to have a good audio quality out of my mixer. But this will be determined later, when I will have all the audio capacities of my audio processor.

You can give your answer in floating point or in integer, because I already know way about how to convert them regarding the need I would have.

I know that this question has been asked many times, but I don't feel like I had the satisfactory answer that I'm seeking. I need a general formula to my specific problem. If I'm finding good answer, I would be very happy to share my simple low level mixer code into a github and link it here.

Edit :
For those that might fall here while searching for an audio mixer. The formula I gave is one of the correct one. To mix samples you'll just have to sum them. Easy as that. Then the clipping part is up to you, and what you want to do : most of the mixer let the user make the sound clip, it allows you more flexibility on the amplitude of each sample, ie : you want the kick to sound 2 times louder than the snare: go for it and find a way so that it won't clip.


texts: Would the ratio of Autocorrelation value and autocorrelation of noise for different data types be different? in Matlab

texts:
Considering a model of the form:

$y[n] = x[n] + noise[n]$ and

$y[n] = u[n] + noise[n]$

where the term noise represents Additive White Gaussian Noise of variance $\sigma^2_{noise}$.
I want to understand if there is any difference in the value of the terms: $\sigma^2_x/\sigma^2_{noise}$ and $\sigma^2_u/\sigma^2_{noise}$ by varying the noise level.

$x$ represents a signal generated from Gaussian distribution whose variance is $\sigma^2_x$. In Matlab, I have used the command rand() to generate $x$.

$u$ represents another signal taking values zero or 1 having variance $\sigma^2_u$. In Matlab, I have used rand()>0.5 to generate thie binary valued data.

Based on my understanding the numerator in the ratio is the autocorrelation of the data mathematically expressed as $E[\mathbf{x}[n]\mathbf{x}[n]^T]$ and  $E[\mathbf{u}[n]\mathbf{u}[n]^T]$. Why are these terms giving same result?

In most papers, the variance of the data irrespective of whether it takes numeric or symbols is considered to be 1.

Question1) Is my approach of implementing the ratio correct? I cannot understand what the value of $\sigma^2_x$ and $\sigma^2_u$ should be in practice. Implementing this way would give the same value for the ratio for both the different kinds of data $x$ and $u$. Would the auto correlation for numeric values x and binary data u be the same value?

Question2) If the ratio is the same then what is the need for using binary signal?

N = 128;
var_x =1;
var_u = 1;
signal_x = rand(1,N);
signal_u = rand(1,N)>=0.5; % this creates 0/1 data
noise = 0;
index = 1;
for noise = 0:5:30
noise_var  = 10^(-noise/10);
ratio_x(index)= 10*log10(var_x/noise_var); %in dB
ratio_u(index)= 10*log10(var_u/noise_var);
index = index+1;
end


texts: What is the meaning of these Fixed Point Properties e.g.CustomSectionInputDataType?

texts:
I was trying tom implement a Butterworth filter and i came across a example on Mathworks site which emphasizes to customize various terms related to the filter.
Following are some of the terms i could not understand and i was also not able to find suitable  info about them in matlab documentation-

1.CustomScaleValuesDataType
2.SectionInputDataType
3.CustomSectionInputDataType
4.SectionOutputDataType
5.CustomSectionOutputDataType

Can anyone please help me to understand this?

This was the link i was referring to-
https://in.mathworks.com/help/hdlfilter/examples/hdl-butterworth-filter.html


texts: Why Fourier Transform doesn't show obvious periodicity?

texts:
I have a clearly 4Hz periodic signal (see image below). However when I apply a fourier transform, there is no peak at 4Hz. Also, when I do a power spectrum density (using Matlab's pwelch algorithm), there is no peak at 4Hz.

My original signal is sampled at 5ms, and is 141 bins (705ms) long.

Why is this & which alogithm will detect the 4Hz peak?

Thanks


texts: FFT of a $N$-length real sequence via FFT of a N/2-length complex sequence

texts:
I would like to optimise FFT algorithm by only using a FFT of size $N/2$ to transform a real sequence of length $N$. There are already many questions and many excellent answers on the matter (here for example).

I've been trying to follow the linked answer to implement the splitting algorithm in Matlab and I checked against the standard FFT Matlab routine. I found a shift of 1 sample in the spectrum (sine input) and a weird peak at a (seems to me) random frequency. I would really appreciate a review of my code or simply a hint to debug the problem.

As a side note, I'm trying this route because I need a memory efficient FFT for an embedded device.

close all;
clear all;
clc;

%% Parameters
sample_rate = 50; % Hz
frequency = 10; % Hz
N = 256;
p = 1:N;
t = (p-1) / sample_rate;

f = sample_rate*(0:(N/2 -1))/N; % frequency vector

%% Generate real sequence
data = sin(2*pi*p*frequency/sample_rate);

figure();
plot(t, data);
title('Input');

%% Split even-odd
% Even indeces are odd because
% Matlab starts indexing from 1
idx_even = 1:2:N;
idx_odd = idx_even+1;
data_even = data(idx_even);
data_odd = data(idx_odd);
data_cpx = complex(data_even, data_odd);

%% Perform FFT of length N/2
X = fft(data_cpx, N/2);

%% Perform Butterfly
k = 0:(N/2 -1);
W = exp(-2*pi*1i*k/N);

X_even =     0.5*(X(k+1) + conj(X(N/2 - k)));
X_odd  = -1i*0.5*(X(k+1) - conj(X(N/2 - k)));
G = X_even + W .* X_odd;

%% Perform FFT of length N to check results
X_N = fft(data, N);

%% Plot
figure();
plot(f, real(X_N(1:128)));
hold on
plot(f, real(G));
title('Real');

figure();
plot(f, imag(X_N(1:128)));
hold on
plot(f, imag(G));
title('Imag');

figure();
plot(f, abs(X_N(1:128)));
hold on
plot(f, abs(G));
title('Mag');

figure();
plot(f, angle(X_N(1:128)));
hold on
plot(f, angle(G));
title('Phase');


texts: how could window function really suppress those frequencies that we don't want?

texts:
I know we can use window function to decrease leakage phenomenon by making the end of a non-periodic signal to be zero at both end, so that to make it like "continuous and periodic" when doing fft.

However, when we multiply window function with original signal
(ex: w(t)xs(t)), it's like doing convolution (w(f)*s(f)) in frequency domain, just like w(f) function slide from left to right and integrate with the area with s(f).

And my question is that since it would still integrate the area which is not the real frequencies of the original signal, and have values after convolution, how could this method really suppress those frequencies?

What make me confused is because I think using window function is the same as the method of doing FIR, which also multiple a window function with the ideal low pass signal in frequency domain. And the FIR characteristics are related to how many areas during integration (ex: the size of transition band).

Thank you very much!


texts: How do I turn my sensor samples into a cleaner, more useful curve?

texts:
This is for one of my fish tanks, which has dimmable LED lights. It's a bit of an experiment where I'd like to see what'll happen if the "sun" above the tank "rises" and "sets" in relation to the actual sun, as in, that glowing ball of fire in the sky. To that end I've hooked them up to a microcontroller reading a light sensor. What works really well is when there's a cloud moving in front of the sun, or the opposite. What's not working so well is to follow the course of day. I've logged the data for a few days. Here's two examples:

(mostly overcast day)

(this one had more sun in it, but also rain showers, and me working in that room during the late evening)
While I'm very happy with the sensor data as such, my various attempts of cleaning it up (taking the average over 10 seconds, 60 seconds, and even up to 3600 seconds) are not really leading to a nice brightness progression of the fish tank LEDs' brightness setting¹. What I'm looking to get in the end is something more like this:

(https://www.greenquarter.biz/green-engineering/daylight-harvesting/)
Is there some mathematical way of working out the general trend of a curve which you know has the general property of first rising from 0 for a while and then falling back to 0 for a similar while, without having to be able to look into the future?²
Quick edit before going to bed:

(Running 60-minute average)

¹No, it has nothing to do with the human eye's lightness constancy, I'm already compensating for that.
²If it matters, the microcontroller I'm using is an Atmel SAMD21, and while it doesn't have an FPU, it's sufficiently fast that I don't mind working with floating-point numbers.


texts: FFT of SIN waves with different phase delays

texts:
I have come across a peculiarity of FFTs which has got me somewhat baffled.
I've simply summed up 101 sine waves and taken the FFT using this matlab script :

clear all
f=1e9;                                 % Centre Frequency 1GHz
df=2.5e6;                                % Carrier Frequency 2.5MHz
Time=linspace(-100e-9,100e-9,1000);      % Region of time to simulate over
delay=0;
Voltage=Time.*0;                        % Initialise Voltages to zero
for loop=-50:50                         % Sum 101 carrier Frequencies
    Voltage=Voltage+sin(2.*pi().*(f+df.*loop).*(Time-delay));
end
figure(1)                               %Plot Time dependent response
subplot(2,1,1)
plot(Time,Voltage)
subplot(2,1,2)                               %Plot Frequency Content
dt=Time(2)-Time(1);
frequency=linspace(-0.5/dt,0.5/dt,1000);
spectrum=fftshift(fft(Voltage));
plot(frequency,abs(spectrum))

The output is as I had expected, with the correct frequency content :

However, if I simply add a significant time delay (by re-running the script with delay=150e-9;, such that the main constructive interference lobe disappears outside the calculation window) the frequency content of the resulting time trace collapses to two peaks.

However the time trace is still the summation of 101 sin waves albeit now out of phase because of the introduced delay ?? Intuitively I would have expected the absolute frequency content of the trace to be preserved and only the phases modified by the delay. Upon reflection I can perhaps understand that the frequency content must be modified on energy conservation grounds, but can anybody rationalise what is going on here ?


texts: Get phase information from bursty time series with amplitude variability and sharp edges

texts:
I have time series like this one, and I am trying to get the phase information from it. If it matters why, because I need to compute the precision of the cycles and amplitude variability is introducing a lot of "noise" to the signal that I want to get rid of.

Time series:

Purpose:
At the end of the day I want to make an autocorrelation of the signal to compute the Quality Factor.

Problem:

The problem is that since the signal has a very sharp edge at the beginning of every cycle, any approach I try to do give me a noisy phase information that I can't use.
Computing the Fourier Transform of the original series also contains the amplitude variability, which I don't want to use.

I tired to approaches, using a Wavelet transform and a Hilbert transform. From any of these to approaches I get a Complex time series from which I can compute the phase, in principle with out amplitude variability.

The problem is that because of the sharp edges of the times series both approaches are failing to get a smooth phase profile of the signal.

Hilbert Transform:

I am using the Hilbert transform of Matlab, which uses FFT. The thing is that since the edges are sharp, FFT is failing the get the correct profile of the signal.
I have the same problem if I make a FFT and then the inverse transform, I get this type of profiles, in which there are maximums where there are not supposed to be, and of course the phase profile, or the cosine, it's not smooth at all.

Wavelet:

Using a Complex Gaussian Kernel, "cgau1" that is the one that work best.

In the figure: Wavelet return and some cuts of the angles.

In the figure: making some cuts for different scales of the angle I get the followings cosine of the phase profile. Which again they are not smooth or they loose cycles, depending on the scale.

Question

How can I get a smooth phase profile, or the cosine of the phase, from a time series with high amplitude variability and with sharp edges?

python or matlab approaches are useful

Thank you very much.


texts: Accurately measuring relative distance between a set of fiducials (Augmented reality application)

texts:
Let's say I have a set of 5 markers. I am trying to find the relative distances between each marker using an augmented reality framework such as ARToolkit. In my camera feed thee first 20 frames show me the first 2 markers only so I can work out the transformation between the 2 markers. The second 20 frames show me the 2nd and 3rd markers only and so on. The last 20 frames show me the 5th and 1st markers. I want to build up a 3D map of the marker positions of all 5 markers.

My question is, knowing that there will be inaccuracies with the distances due to low quality of the video feed, how do I minimise the inaccuracies given all the information I have gathered?

My naive approach would be to use the first marker as a base point, from the first 20 frames take the mean of the transformations and place the 2nd marker and so forth for the 3rd and 4th. For the 5th marker place it inbetween the 4th and 1st by placing it in the middle of the mean of the transformations between the 5th and 1st and the 4th and 5th. This approach I feel has a bias towards the first marker placement though and doesn't take into account the camera seeing more than 2 markers per frame.

Ultimately I want my system to be able to work out the map of x number of markers. In any given frame up to x markers can appear and there are non-systemic errors due to the image quality.

Any help regarding the correct approach to this problem would be greatly appreciated.


texts: Faster Algorithm to convolve/correlate two sparse 1-D signals in python (or any language)

texts:
I have two signals which I need to correlate or convolve. Each signal is sampled non-uniformly and the values of the signal I have with me are the timestamp and the magnitude of the signal at that timestamp. The values of the signal at all the other times can be assumed to be zero. The timestamps of the signal have a resolution running in microseconds.

An example of how the signal looks is shown below:

As can be seen, the resolution of the signal is in microseconds and signal is mostly sparse.

If I were to the convolve two signals of this type, I would first have to pad the signals with zeros (since I would have to discretise the signal). While the padding can be done with resolution of microseconds, The number of values to be multiplied becomes too big and the operation becomes increasingly slow. Most of the multiplications in this convolution would be multiplication of zeros(which are pretty much useless). I have therefore chosen a round off value of 2 places (0.xxxxxx becomes 0.xx),since I have to perform 40,000 similar convolutions. I have written my resampling function as shown below.

import numpy as np
import math

def resampled_signal_lists_with_zeros(signal_dict, xlimits):
  '''
      resamples the given signal with precision determined by the round function.
      signal_dict is a dictionary with timestamp as key and signal magnitude is the value of the key.
      xlimits is an array containing the start and stop time of the signal.
  '''
  t_stamps_list = list(signal_dict.keys())
  t_list = list(np.arange(int(math.floor(xlimits[0])), int(math.ceil(xlimits[1])), 0.005))
  t_list = [round(t, 2) for t in t_list]
  s_list = list()
  time_keys = [round(t, 2) for t in t_stamps_list]

  i = 0
  for t in t_list:
    if i < len(t_stamps_list):
      if t==time_keys[i]:
        s_list.append(signal_dict[t_stamps_list[i]])
        i+=1

      else:
        s_list.append(0)
    else:
      s_list.append(0)
  return t_list, s_list

The correlation of two signals padded in the above manner is done using scipy as follows:

from scipy.signal import correlate
output = correlate(s_1, s_2, mode='same')

The output calculated in the above manner is pretty slow .Since the signal is pretty sparse and most of the multiplications in the signal are multiplications of zero, I think there should be a better way to do the same operations. Is there a way to get the result of the convolutions of the two sparse signals faster?


texts: Matlab group delay not corresponding to reality

texts:
I studying effects of digital filters group delay (IIR). In Matlab there is a function grpdelay which allows for particular frequency get a value of delay in samples
But after filter applying and comparing of pictures there are diffent results. In rude words if I expect 50 samples delay under grpdelay analysis, in fact I see 100 samples, If I expect 30, I see 70 and so on. What I doing wrong? I use simple sin processed gauss window for example
clear;
close all;

t = 0:0.008:600;

x = sin( 2 * pi * t * 1.3 ) .* gausswin(length(t) )';
x = [ zeros(1, 1000) x zeros(1, 1000) ];

fs = 125;
fc_low = 0.5;
fc_high = 3;
order = 4;

[b,a] = butter(order, [ fc_low/(fs/2), fc_high/(fs/2) ] );
y = filter(b, a, x);

plot(x)

hold on;

plot(y)

legend('orig', 'filtered')

figure;

grpdelay(b, a, 512, fs);

below there is example of theoretical and model misfitting

So, the result is that Matlab function grpdelay gives not accurate results or I have errors in calculation/interpretation ?


texts: Ratio of expected values of squares of errors in power quantities, in dB

texts:
An amplitude quantity $a$ is estimated by two estimators $\hat a_1$ and $\hat a_2.$ The error of estimator $\hat a_2$ is compared to that of $\hat a_1$. One possible comparison is the ratio of expected values of absolute error:

$$\frac{\operatorname{E}[|a - \hat a_2|]}{\operatorname{E}[|a - \hat a_1|]}.$$

This is a ratio of amplitude quantitities and can be expressed in decibels as:

$$20 \log_{10}\left(\frac{\operatorname{E}[|a - \hat a_2|]}{\operatorname{E}[|a - \hat a_1|]}\right) \text{ dB.}$$

Another possible comparison is the ratio of expected values of square error:

$$\frac{\operatorname{E}[(a - \hat a_2)^2]}{\operatorname{E}[(a - \hat a_1)^2]}.$$

This is a ratio of power quantities and can be expressed in decibels as:

$$10 \log_{10}\left(\frac{\operatorname{E}[(a - \hat a_2)^2]}{\operatorname{E}[(a - \hat a_1)^2]}\right) \text{ dB.}$$

The squares of the estimators can be used to estimate $a^2.$ The estimators $\hat a_1^2$ and $\hat a_2^2$ can be compared by the ratio of expected values of square error:

$$\frac{\operatorname{E}[(a^2 - \hat a_2^2)^2]}{\operatorname{E}[(a^2 - \hat a_1^2)^2]}.$$

This is a ratio of squares of power quantities. How should it be expressed in decibels?

$$10 \log_{10}\left(\frac{\operatorname{E}[(a^2 - \hat a_2^2)^2]}{\operatorname{E}[(a^2 - \hat a_1^2)^2]}\right) \text{ dB,}$$

or perhaps:

$$5 \log_{10}\left(\frac{\operatorname{E}[(a^2 - \hat a_2^2)^2]}{\operatorname{E}[(a^2 - \hat a_1^2)^2]}\right) \text{ dB} = 10 \log_{10}\left(\frac{\sqrt{\operatorname{E}[(a^2 - \hat a_2^2)^2]}}{\sqrt{\operatorname{E}[(a^2 - \hat a_1^2)^2]}}\right) \text{ dB},$$

where the square roots are power quantities, or is there no commonly agreed convention?


texts: Why does a simple delay system with $h[n]=\delta[n-n_0]$ have a linear phase?

texts:
Let's consider a simple delay system with impulse response $$h[n]=\delta[n-n_0]$$ We also know that it has linear phase with group delay of $n_0$
.

In general, a FIR filter with symmetric (or anti-symmetric) coefficients with filter order $M$ always has linear phase and its (group) delay is $M/2$. Now my question is:
1) What is the filter order of the simple delay system?
    My guess is that the order is $n_0$ (since it is number of unit-delays)
2) Does it have symmetric filter coefficients?
    Apparently it is not, but not quite sure.
3) Why does it have a linear phase?


texts:
When a digital device (PC, portable media player, etc.) plays a digital audio file (ogg, mp3, flac, etc.) the audio signal output is always the same, regardless of the device type/brand, am I right? Thus the sound quality should not differ. For instance, I have an iPod that plays an mp3 file. If I play the same file on a different brand media player the sound quality should be the same since it is the very same signal (no equalizers, or any sound alterations).

Is this true? If it is true then I reckon it is only the headphones/speakers that matter as far as the sound quality is concerned.

P.S: all of the answers are really nice! wish i could accept 'em all


texts: Time-Delay beamforming delay computation

texts:
I'm trying to implement a beamforming algorithm in python and I'm stuck with the time-delay computation for a rectangular array.

I'm using the following equation to compute the delays based on the desired steering angle :
$$d_{i,j} = \frac{1}{c}\sin(\theta)(\cos(\phi)x_{i,j} + \sin(\phi)y_{i,j})$$

Where $d_{i,j}$ is the delay for the element $(i,j)$, $(\theta,\phi)$ is the azimuth/elevation steering angle, and $x_{i,j},y_{i,j}$ are the coordinates of the element $(i,j)$.

I am using Matlab's PhasedArray toolbox as a ground truth to verify my work, but I'm getting different results.

Here is the Matlab code I use to get the time delays :

array = phased.URA('Size', [4,4],...
                   'ElementSpacing', 0.042);
delay = phased.ElementDelay('SensorArray',array, 'PropagationSpeed', 340.0);
tau = delay([-30;-20])*fs

And here the Python equivalent :

ang_dft = np.array([-30, -20], float)

myArray = pa.URA(4, 4, 0.042, 0.042)
TDBeamformer = bf.TimeDelayBeamformer(myArray, ang_dft, 340.0, 8000.0)

myArray.plot_array()
TDBeamformer.compute_delays()

print(TDBeamformer.sec2spl(TDBeamformer.delay))

plt.show()

Where the compute_delays method is the following :

def compute_delays(self):
    for x in range(0, self.array.sizeX):
        for y in range(0, self.array.sizeY):
            self.delay[y, x] = (1/self.c *
                                np.sin(np.pi / 180. * self.steeringAngle[0]) *
                                (np.cos(np.pi / 180. * self.steeringAngle[1]) * self.array.coordinatesX[y, x] +
                                 np.sin(np.pi / 180. * self.steeringAngle[1]) * self.array.coordinatesY[y, x]))

Could you help me pinpoint where the problem lies? Is my equation wrong (I found it from this website)? Is it a geometry problem in the code?


texts: Find frequency/phase response of a microphone

texts:
I am trying to find the frequency/phase response of two cheap microphones bought off Ebay. Basically, I want to graph the frequency responses and phase responses of both microphones on the same axis so I can compare their sensitivity in MATLAB.

One approach I have read is downloading and playing a file that consists of pink noise through a speaker and simultaneously recording this using both microphones, then using the two recorded signals and the original pink noise file to get the bode plots. This is achieved by first taking the FFT of the recorded signals and of the original pink noise signal, and then using the bode function in MATLAB. The frequency/phase response of the speaker is unknown, however I'm more interested in the differences between the two microphones, so my thoughts are this wont affect my analysis. Any thoughts on this?

Is there an easier approach to finding the frequency/phase response of two microphones that I am not aware of using MATLAB?


texts: Detecting blobs that connects to any other blob, maybe with OpenCV

texts:
In the image, I connects to universe, but II and III not.

I need to detect both II and III, also I if possible.

Is it possible with current computer vision libraries?


texts: Filter Banks and FFT, are they similar conceptually?

texts:
I'm using Filter Banks (Narrow Band Pass) created using the functions in Matlab. However,  I wanted to know if there are fundamental limitations like that of FFT which also affect the Filter Banks. Is there any co-relation between the two?

Basically, I don't know if they share the same problem of time and frequency being inversely related. If they do share this problem, this means that we cannot have higher time and higher frequency resolution at the same time.


texts: How much can we downsample by? and how much can we upsample by?

texts:
e.g. let's say we have a signal with B = 100 Hz that is sampled at 400 Hz. The anti-aliasing filter will remove replicas, and there is nothing to alias, and we downsample. If we downsample by 2, we will get signal of 200 Hz. But what's stopping us? Can we downsample by 4, or 8?

(why can't we downsample by 4, get B = 400Hz, then have a new sampling rate at 800 Hz?)

Similarly, is there a limit to how much we want to upsample?


texts: FIRFilter with low cut-off compared to sampling frequency

texts:
Sorry if my question is some FIRFiltering basics, but I couldn't find a specific answer after some search.

I am trying to design a decimating low pass filter in C (actually in OpenCL with GPU) with the FIR approach by doing convolution.
I have a software that generates the coefficients based on the cut-off frequency and the sampling frequency.

I realized that my filter generally works if the cut-off frequency is not too far away from sampling frequency, for example fsample = 1e7 and fcutoff = 1e4 would work. (maybe the attenuation is not enough). But if I increase the fsample to 1e8, the attenuation will be negligible.

Reading from this link:
How many taps does an FIR filter need?
I understand that I will need a number of tap of order 1e5 or 1e6 to have sufficient attenuation (~10 - 20 dB) for my application.

So my question is, in general, how should I design my filter to meet my requirement, if I really have to cut-off something massively oversampled?
If not FIR, could you point me any other better approach, for example, can IIR work better?

Maybe I can follow the approach of multi-stage decimating FIR filter in this page? Great if you provide me with some references regarding how the stages should be determined.
FIR-Decimation and Low-pass filter (taps vs number of input points vs number of decimation stages)

Thanks.


texts: Inverse fourier transform of Hermitian function, getting an imaginary part

texts:
in the cartoon below

It shows that if we take the inverse Fourier transform of a Hermitian function, real part even and imaginary part is odd we should get a purely real function in the time domain.

I tried to replicate this by taking a frequency response I have, zero padding it making it even for the real and odd in the imaginary. Shown below

Although my ifft in Matlab I still have a significant imaginary component, the real is shown left and imaginary on the right.

I also wrote my own script for DFT/IDFT and tested it on a sine function to make sure I get the same as Matlab. Maybe I have missed something?


texts: Destructive interference on the Autocorrelation of the time signal of periodic CPFSK-signals.

texts:
An FSK-signal with a frequency-shift $\Delta F$, a symbol-stream $x$, $x(n) \in {\{-1,1}\} $ and symbol duration $T$ has the complex envelope:

$$f[x](t)=\exp\left(2 \pi i \Delta F T \left(\sum_{l=0}^{j-1} x(l)\right)+\left(\frac{t}{T}-j\right)x(j)\right)\quad\text{for}\quad jT\leq t\leq (j+1)T$$

I am examining periodical symbol-streams for FSK-signals and I found a peculiar dependency between the strength of the periodicity and the modulation index,which I am very unsure of:

The phase of the samples at the start of a symbol is

$$2\pi i \Delta F T \left(\sum_{l=0}^{j-1} x(l)\right)=:\varphi(j,T).$$

So the phase deviation between the samples at the beginning of two consecutive symbols is given by

$$|(\varphi(j,T) -\varphi((j-1),T)|=\pi \eta, \tag{1}$$

where $\eta=2 \Delta F T$ is the modulation index of the FSK signal.

Due to $(1)$ for integer modulation indices $\eta$ the phase at the beginning of a symbol stays constant. Accordingly, a complex autocorrelation of a periodic signal shows a peak at right at the length of a period of the signal.

But, if the modulation index is not integer, the starting phase of the symbols shifts throughout the signal, i.e. $2\pi \nmid \pi \eta $.
A complex autocorrelation of same period signals, shifted in phase, does not show the same peaks and therefore does not exhibit periodicity in the same magnitude.

Is my suspicion  really sound for continuous phase frequency shift keying signals?
Does this hold for MSK-signals too ?


texts: Difference in frequency spectrum of continuous time signal and discrete time signal?

texts:
I am studying chapter 4 of book signal processing first

There i came across a concept which says that(what i have understood)the spectrum of continuous time signal has limited number of lines while spectrum of discrete time signal has infinite number of lines

Why?what is the reason for this difference?

https://imgur.com/a/6IENWha

above link
Is snap of  frequency spectrum diagrams of both continous time and discrete time signals

below link is relevant theory of book
https://imgur.com/a/mH7L7D6


texts: Applications for adaptive and non adaptive histogram equalization

texts:
I have a question regarding histogram equalization and adaptive histogram equalization.

In which circumstance would it be preferable to use adaptive histogram equalization over histogram equalization and similarly, when would it be preferable to use histogram equalization?

Which applications do they find more use in compared to one another?

From using Matlab, I understand the histogram equalization provides a uniform pixel region. And that adaptive histogram locally enhances the contrast on a specific pixel region.

Any thoughts? Thanks.


texts:
I'm new signal processing (this is my first pratical application), so I ask sorry if my question is simple.

My challenge is to see the displacement of accelerometer when I have vibration too. I do not know the vibration completely. The accelerometer will be attached to hoop of started vehicle.

I know that my "displacement have very low frequency" (maybe this statement is not adequate) in relationship vibration. The movement will be do moving accelerometer with hand. The vibration is coming of the vehicle engine.

I will use adxl362, it can work up to 400 Hz and can be programmed with low pass filter of 1/2 or 1/4 of output data rate.

I would like help with steps that I should follow to discriminate if  am vibrating or moving.

PS: 1 - I no need precision.
      2 - After earch detect, the vibration and movement will be stopped.

Tnx very much for any help!


texts: Why plotting the fft result of a signal from -f to f gives different result compared to 0 to f?

texts:
When I plot the fft result of an audio signal on the positive side of frequency axis I get the following figure, I understand the signal on the right side of the Nyquist point is just a mirror image of the original signal, the sharp noise is somewhere around 1700Hz:

Now when I plot the fft of the same signal from negative frequency to positive, that noisy component is now sitting somewhere around 2300Hz and when I check the fft of the audio signal, the sharp noise is indeed at 2300Hz and not 1700Hz, the audio content also is no longer below 1KHz and it's now instead just below 4KHz. What exactly is going on here ? Why the top plot appears to be giving me the wrong information ?

I have tried to find the reason online but I don't think I am googling the right words, hence can't find anything.


texts:
Actually my question is how to get the weighted average of two or three IIR filters given their transfer functions, which is nearly equal to this question. How to calculate the total transfer function of the summation of two IIR filters? For example, I have two filters as following
$$H_1(z) = \frac{b_0+b_1z^{-1}+b_2z^{-2}}{1+a_1z^{-1}+a_2z^{-2}}$$
$$H_2(z) = \frac{b_0'+b_1'z^{-1}+b_2'z^{-2}}{1+a_1'z^{-1}+a_2'z^{-2}}$$
I can directly solve the sum of these two fractions as $H(z) = H_1(z)+H_2(z)$, but the problem is that my filter orders are generally around dozens, I have no idea how to do it.
I also know that I can put a signal through this two filters and sum the output together, but it requires double computational complexity.
Is there any other way to solve this problem? Thanks!


texts: Creating a matrix of Gaussian wavelets at dyadic scales

texts:
I need to create a diagonal matrix containing the Fourier coefficients of the Gaussian wavelet function, but I'm unsure of what to do.

Currently, I'm using function ConstructHaarWaveletTransformationMatrix to generate the Haar wavelet matrix and am taking the rows at dyadic scales (2, 4, 8, and 16) as the transform:

M = 256
H = ConstructHaarWaveletTransformationMatrix(M);
fi = conj(dftmtx(M))/M;
H = fi*H;
H = H(4,:);
H = diag(H);

etc.

How do I repeat this for Gaussian wavelets? Is there a built-in MATLAB function which will do this for me?

For reference, I'm implementing the algorithm in section 4 of the paper in Compressed Sensing for Wideband Cognitive Radios.

I need to do edge detection on a 1D signal, which is constructed like so:

edges = [50, 120, 170, 200, 220, 224, 256];
levels = [24,3,30,0,36,0,0];
idxs = zeros(1,max(edges));
idxs(edges(1:end-1)+1) = 1;
psd = levels(cumsum(idxs)+1);

When I just wavdec and the differentiation matrix:

[y1,y2] = wavedec(psd, 0, 'haar');
z0 = gamma*y1.';

I get exactly what I want.

The point, however, is that I want to recreate this edge information compressively (that is, by solving a linear program).

The fact that using the matrix above creates gibberish is my current sticking point. I want to replace it with something else and see if it still works.


texts:
I'm a computer science student and I took a network course this semester which assumes some signal processing backgrounds which I have not. Last time We talked about ASK.

ASK is Amplitude Shift Keying,

it uses two different amplitudes to describes $0$ or $1$. We can use several levels to represents more symbols.

$$A\sin(2\pi fnt+\phi)$$

How many harmonics does a signal with an ASK coding has?

If I have such a message to broadcast: $01101101$

How do we code in ASK this series of bits when the frequency used is: $f=1kHz$? And when the bit-time is $T_b=1ms$?

And what about varying its amplitude? If we have 4 levels for instance?


texts: Fast Recursive 1D Signal Smoothing - IIR / Auto Regressive Implementation of Gaussian Smoothing

texts:
I have just begun to dive into the field of signal processing, but there is the need to program a digital filter, that has to smooth a realtime signal from a sensor device. As far as I know, in my case a FIR filter is not suitable, because I have a lot of samples (1000) per sec and the convolution that is needed for appling a gaussian smoothing takes to long. Therefore I need some IIR Filter, that is capable of deliver results, that are comparable with gaussian filtering. I have already tested the first order lag filter, but this filter does not smooth out the signal if it was constant over time before.

Do you have any suggestions where to start? What are suitable fast recursive smoothing filters that look like gaussian smoothing? Some hints will be usefull for my to get the keywords for further reading.


texts: What is the interpretation of Fourier Transform containing only imaginary part?

texts:
The FT of a unit step function is taken as:

$$
X(\omega) = \int_0^\infty e^{-j\omega t}dt = \frac{-1}{jw}e^{-j\omega t} \Biggr |_{0}^{\infty} = \frac{j}{\omega}
$$

The transform only has the imaginary component and nearly (actually) all literature that I have read deem this transform as non-existent because of the missing real component.

Unfortunately, I have not found an intuitive explanation that describes why this transform is invalid/non-existent.


texts: About convergence of KL divergence: if the two probability distributions are type, does the law of large number work?

texts:
If I pick $N$ samples from $P_X$ and $P_Y$, they are two independent discrete distributions.

$X_1,X_2,\ldots,X_N$ are drawn i.i.d from $P_X$, and $Y_1,\ldots,Y_N$ are drawn i.i.d from $P_Y$.

I got $P_X'$ is the type of $X_1,X_2,\ldots,X_N$  and $P_Y'$ is the type of $Y_1,Y_2,\ldots,Y_N$.
That is, $P_X'$ converges to $P_X$ almost surely as $N$ goes to infinity, and
$P_Y'$ converges to $P_Y$ almost surely as $N$ goes to infinity.

My claim is $D(P_X'\|P_Y')$ will converge to $D(P_X\|P_Y)$ as $N$ goes to infinity too, in the above intuition.

But as the KL-divergence $D(\cdot\|\cdot)$ is not continuous, how to start a rigorous proof that the claim is true or false?

Because there are something in log in KL-divergence, It's hard to write down a form like strong law of large number in wikipedia.


texts: Weird results when using MATLAB's $\tt butter$ function?

texts:
I'm trying to learn some of the basic MATLAB filter design functions and have been experimenting with the butter function. This function calculates the transfer function coefficients of a Butterworth filter given a desired filter order and normalized cutoff frequency.

I have previous experiences with designing prototype filters and then modifying them, however, I'm having a difficult time understanding MATLAB's definition of normalized frequency.

I'm used to designing Butterworth filters with a normalized frequency of $1 \textrm{ rad/s}$, yet MATLAB insists that the normalized frequency has something to do with a sample rate? I assumed the butter function was used for designing analog filters?

Furthermore, when I plot the frequency response of the function, the cutoff frequency I specify (Wn parameter) never seems to correspond to the $-3\textrm{ dB}$ mark in the plots.

Any thoughts on these?


texts: Implementation of goertzel algorithm in C

texts:
I am implementing BFSK frequency hopping communication system on a DSP processor. It was suggested by some of the forum members to use Goertzel algorithm for the demodulation of frequency hopping at specific frequncies. I have tried implementing  the goertzel alorithm in C. the code is follows :

float goertzel(int numSamples,int TARGET_FREQUENCY,int SAMPLING_RATE, float* data)
{
int k,i;
float   floatnumSamples;
float   omega,sine,cosine,coeff,q0,q1,q2,result,real,imag;

floatnumSamples = (float) numSamples;
k = (int) (0.5 + ((floatnumSamples * TARGET_FREQUENCY) / SAMPLING_RATE));
omega = (2.0 * M_PI * k) / floatnumSamples;
sine = sin(omega);
cosine = cos(omega);
coeff = 2.0 * cosine;
q0=0;
q1=0;
q2=0;

for(i=0; i<numSamples; i++)
{
    q0 = coeff * q1 - q2 + data[i];
    q2 = q1;
    q1 = q0;
}
real = (q1 - q2 * cosine);
imag = (q2 * sine);
result = sqrtf(real*real + imag*imag);
return result;
}

When I use the function to calculate the result at specific frequencies for a given dataset, I am not getting the correct results. However, if I use the same dataset and calculate the goertzel result using MATLAB, then I get the results perfectly. I am implemented the algorithm using C, with the help of some online tutorials that I found over the internet. I just want to get the view of you guys if the function is implmenting the goertzel algorithm correctly.

Thanks,
Anshu


texts: Strange behavior from frequency to angular frequency in FFT

texts:
I have written a code to compute the Fast Fourier Transform of a simple complex exponential with frequency $f=50.0$, using scipy.fft. The code is written below:
import numpy as np
from scipy.fft import fft, fftfreq, fftshift
import matplotlib.pyplot as plt

T = 10
dt = 0.001
f = 50.0

x = np.arange(0.0, T, dt)
y = np.exp(1.j * 2 * np.pi * f * x)

xf = fftfreq(len(x), dt)
xf = fftshift(xf)

yf = fft(y)
yf = fftshift(yf)

plt.plot(xf, 1.0/len(x) * yf, label='yf')
plt.plot(xf, 1.0/len(x) * np.abs(yf), label='abs(yf)')
plt.xlabel('f')
plt.legend()
plt.show()

The output of the calculation is a frequency spectrum with a single peak, for both yf and its absolute value (which are perfectly overlapped) centered at frequency $f=50.0$:

Now instead of using the frequency, I assume to work with an angular frequency $\omega=50.0$:
import numpy as np
from scipy.fft import fft, fftfreq, fftshift
import matplotlib.pyplot as plt

T = 10
dt = 0.001
omega = 50.0

x = np.arange(0.0, T, dt)
y = np.exp(1.j * omega * x)

xf = fftfreq(len(x), dt)
xf = fftshift(xf)

yf = fft(y)
yf = fftshift(yf)

plt.plot(2 * np.pi * xf, 1.0/len(x) * yf)
plt.plot(2 * np.pi * xf, 1.0/len(x) * np.abs(yf))
plt.show()

In the plot I have multiplied xf by $2 \pi$ to convert from frequency $f$ to angular frequency $\omega$.
Surely this corresponds to a frequency $f = \frac{\omega}{2 \pi}$ different from the first example but I would have expected similar behavior. In this case yt have a kind of dispersive shape, while its absolute value is very broad, it is not normalized and it is not perfectly peaked at $\omega=50.0$ but it is slightly shifted to $50.26$.

Could someone give me some explanation about this unexpected behavior?
If I have to work with raw data, coming from a complex calculation, expressed in terms of angular frequency $\omega$, how can I obtain trusty results from the second programs?


texts: Looking for algorithm to extract line network from binary images within LabView

texts:
I would like to create a LabView VI that receives an input image and the RGB-color codes of the powerlines depicted in the image. The output should (eventually) consist of id'd polygonial connections and their coordinates.

Right now i've finished creating a VI that converts the source image to a binary image with the respective line color types extracted into a binary image (see below).

I've been trying to apply some of edge/line finding VIs supplied by LabView Vision, but the results have not been very convincing.

I have considered implementing some sort of vectorization algorithm like 'Straight Line Extraction Using Iterative Total Least Squares Methods' or even some sort of line-walking-and-mapping algorithm, but before I do so I'd like to know if there already is some kind of fitting algorithm out there.

Attached you can find examples for the source and binary images.

Thanks in advance for your input!


texts: Can I use FFT to interpret accelerometer gestures?

texts:
I have the need to detect two different gestures that occur when looking at accelerometer data. Here's a run down (as brief as I can make it):

Lets say an iPhone is being oscillated back and forth while remaining face-up. The  user can either do one oscillation (once forward or backward, gesture one) or a continuous oscillation for any amount of time (gesture two).

Currently my code maintains a list of the device's accelerometer data (y-axis) from the last 50 frames . Each frame this data is put through an FFT algorithm (this one http://goo.gl/yi3mn), and then I try to interpret the frequency domain given. I have noticed a strong relationship between the speed of the oscillations, and the power of the mid to low range of the frequency space.

The problem is that I need to detect (in real-time as the movement animates something on screen) whether the just finished oscillation is a singular one, or one that continues into another oscillation in the opposite direction. Here I should point out that accelerometers work with pressure plates. When an oscillation stops the input data will show the opposite value from when the oscillation was being performed. This makes it difficult to see any obvious difference between a single forwards movement and a double forwards then backwards movement (at the end of the each first oscillation).

Can anyone suggest how I might make use of an FFT to decipher an accurate magnitude value (or something else useful) that I can use to tell the difference between the gestures described above?

I have recorded the data I am using to text files. It is the raw data (g-force values) along the y-axis of the accelerometer input.

http://pastebin.ca/2108123 shows data for 2 singular oscillations (I hold the device still at the start, end, and between the two oscillations).

NOTE: new raw data with 20 oscillations has been uploaded, but not yet plotted.
http://pastebin.ca/2108387 shows data for 20 continuous oscillations (I hold the device still at the start and end).


texts: Derivation of the lowpass to bandpass transformation

texts:
I have a basic question.
The "well known" lowpass to bandpass transformation is
$$
s \longmapsto \frac{\bar{s}^2 + \omega_1\omega_2}{\bar{s}(\omega_1 - \omega_2)},
$$
which gives a bandpass transfer function of
$$
\frac{1}{s + 1} \longmapsto \frac{\bar{s}(\omega_1 - \omega_2)}{\bar{s}^2 + \bar{s}(\omega_1 - \omega_2) + \omega_1 \omega_2}.
$$
My intuition is that a bandpass should be the product of a lowpass and a highpass.  However, this product gives a different transfer function:
$$
\frac{\omega_1}{s + \omega_1} \frac{s}{s + \omega_2} = \frac{\omega_1 s}{s^2 (\omega_1 + \omega_2) s + \omega_1 \omega_2},
$$
which indicates that the bandpass transformation does not give this cascade of lowpass and highpass.

My question is, how is the bandpass transformation designed, in terms of either combining lowpass filters or by pole placement?

Related question, but using a different derivation technique, and reference is made to the lowpass/highpass derivation, but it is not shown: How is the lowpass to bandpass transformation derived?


texts: Multitapering in time domain to obtain smoother filtered signals?

texts:
Suppose we want to approximate the instantaneous power $z(n) = |y(n)|^2$ of the discrete-time signal $y(n)$, where $y(n)$ is the result of filtering $x(n)$ with a given window $h(n)$, $y(n) = x(n) \ast h(n)$, so $z(n) = |x(n) \ast h(n)|^2$. We need to reduce the fluctuation observed in $z(n)$ (i.e. we want to have a smooth $z(n)$), and to achieve this goal we are given the choice of choosing the window function $h(n)$ to convolve with $x(n)$.

Traditional choices for window functions in signal processing are Gaussian, Hanning, Hamming, Kaiser, Triangular, the list goes on. In applications of spectrogram estimation, it is also common to employ the Multitaper approach to reduce the variance of the estimated spectrum, i.e., by averaging out the results of computing $k=1,...,K$ spectrograms, the obtained mean spectrum is smooth.

Could we use the multitaper technique in time domain to obtain smooth approximations of $z(n)$? More specifically, if $h_1(n),...,h_K(n)$ are $K$ orthonormal tapers, shouldn't

$z(n) = \displaystyle \frac{1}{K}\sum_{k=1}^{K}|x(n) \ast h_k(n)|^2$

act like an average filter, yielding to a smoother $z(n)$ than what would be obtained by using a sole window $h(n)$, i.e. $z(n) = |x(n) \ast h(n)|^2$?


texts: translating 1 dB increase in thermal noise to interference to noise ratio (I/N)

texts:
I often see statements like "1 dB increase in the receiver thermal noise density results in an interference-to-noise density ratio of -6 dB". How can I obtain this result mathematically? Assuming the noise power is $N$ (linear domain) and $N'$ (after the 1 dB increase to $N[\text{dB}]$), I have written the following. But I can't get to a -6 dB term?
\begin{align}
10\text{log}\bigg(\frac{I}{N'}\bigg) &= 10 \text{log}(I)-10 \text{log}(N') \\
&= 10 \text{log}(I)-N'[\text{dB}] \\
&= 10 \text{log}(I)-(N[\text{dB}]+1) \\
&= 10 \text{log}(I)-(10 \text{log}(N)+1) \\
&=10 \text{log}\bigg(\frac{I}{N}\bigg)-1
\end{align}


texts: Deriving the Matrix Inversion Lemma for RLS Equations vs the Woodbury Derivation

texts:
Can any one help me in deriving the matrix inversion lemma rule for RLS  algorithm?

I don't know how to start with. Many books have just stated but they haven't derived it.


texts: Set of all images with a certain MSE between themselves and a given image

texts:
Mean squared error between two $M\times N$ images, $I[m, n]$ and $D[m, n]$ is defined as:
$$
MSE(I, D) = \frac{1}{M\times N}\displaystyle \sum_{m=1}^{M}\sum_{n=1}^{N}(I[m, n]-D[m, n])^2
$$
For a given image, $I[m, n]$, and a particular value, $\epsilon $, how to calculate all images, $D[m, n]$, that:
$$
MSE(I, D) = \epsilon
$$


texts:
So I need to make a function that will be called as function [t, pulse] = PAM (array). My function should have a result like this when I use array = [1 0.3 0.8 0.6], which array is the width of my digital pulsing signal.

Also please not a complete answer, I think it's better to think alone, but I'm in a point that I don't know what to do.


texts:
I just have a small theoretical question that buggers me to be honest.

I have a sound file that I found its spectrum on linear and on logarithmic scale. Alright, I have done this, it's a plain standard procedure. The problem is that I cant really understand the plot on logarithmic scale. (it just looks like a typical logarithmic plot with dB, and it is correct by the way) .

Well, the visual differences are obvious, the linear plot frequencies are quite scattered and they peak in a non homogenous way, while I notice that on log scale the frequencies peaks are very close to each other and they don't seem to have so many extreme values. Could someone explain me this a bit ?

Cheers


texts: Can we determine if a filter is butterworth or chebyshev from its physical topology?

texts:
Can all analog filters be classified as Bessel, Elliptic, Butterworth or Chebyshev?

Given a physical ladder topology of several stages of {L, C or LC} in {series or parallel}, is it possible to always classify it to one of the above?

(I got this question while trying to design a steep filter by trial and error on spice, and then realized that it has too long response, as a Chebyshev filter would have, so I wanted to optimize it.)